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Summary 

A  finite-difference  method  is  used  to  investigate  compressible,  laminar 
boundary-layer  flows  of  a  dilute  dusty  gas  over  a  semi-infinite  flat  plate. 
Details  are  given  of  the  implicit  finite-difference  schemes  as  well  as  the 
boundary  conditions,  initial  conditions  and  compatibility  conditions  for 
solving  the  gas-particle  boundary-layer  equations.  The  flow  profiles  for 
both  the  gas  and  particle  phases  were  obtained  numerically  along  the  whole 
length  of  the  plate  from  the  leading  edge  to  far  downstream  of  it.  The 
finite-difference  solutions  in  the  large-slip  region  and  the  small-slip 
region  are  compared  with  the  asymptotic  solutions  and  good  agreement  is 
achieved.  The  boundary-layer  characteristics  of  interest,  including  the 
wall  shear  stress,  the  wall  heat-transfer  rate  and  the  displacement 
thickness,  are  calculated.  The  alteration  of  the  flow  properties  owing  to 
the  presence  of  particles  is  discussed  in  detail.  It  was  found  that  the 
boundary-layer  flow  of  a  dusty  gas  can  be  divided  into  three  distinct  flow 
regimes  which  are  characteri zed  by  quasi -frozen,  nonequilibrium  and 
quasi-equilibrium  flows  and  that  at  a  critical  distance  from  the  leading 
edge  the  particle  velocity  at  the  wall  decelerates  to  zero  and 
near-equil ibrium  is  achieved  between  the  gas  and  particle  flows.  For  the 
laminar  boundary  layer  of  a  dusty  gas,  the  shear  stress  and  the 
heat-transfer  at  the  wall  are  increased  and  the  displacement  thickness  is 
decreased  compared  with  the  pure-gas  case  alone. 
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coefficients  of  the  finite-difference  equations,  Eqs. 

(3.12)-(3.20)  and  ( 3 . 26 ) - ( 3 .27 ) 

grid  parameters  for  the  six-point  difference  scheme,  Eq. 

(3.14) 

coefficient  for  fitting  a  polynomial  to  the  gas  temperature 
near  the  wall,  Eq.  (5.11) 

coefficient  for  fitting  a  polynomial  to  the  gas  velocity 
near  the  wall,  Eq.  (5.10) 


coefficients  of  the  finite-difference  equations,  Eqs. 

(3.12)-(3. 20)  and  (3.26)- (3.27 ) 

grid  parameters  for  the  six-point  difference  scheme,  Eq. 

(3.14) 

coefficient  for  fitting  a  polynomial  to  the  gas  temperature 
near  the  wall,  Eq.  (5.11) 

coefficient  for  fitting  a  polynomial  to  the  gas  velocity 
near  the  wall,  Eq.  (5.10) 
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K  =  Ayn/4yn_i 

heat  conductivity  for  the  gas  phase 
grid  line  in  the  y-direction 

grid  point  at  the  outer  edge  of  the  boundary  layer 
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flow  Reynolds  number  based  on  the  particle  equilibrium 
length,  Re^  =  p*p*>£/nS, 

temperature 

tangential  velocity  in  the  x-direction 

normal  velocity  in  the  y-direction 

function  representing  any  flow  property,  u,  v,  p  or  T 

coordinate  along  the  wall 

coordinate  normal  to  the  wall 

ratio  of  the  specific  heats  for  the  two  phases,  a  =  cjj/c| 
mass  loading  ratio  of  the  particles  in  the  freestream, 

e  =  pf./pi 

displacement  thickness  of  the  boundary  layer 
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density 

density  of  the  particle  material 
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critical  point 

grid  line  in  the  y-direction 
grid  line  in  the  x-di recti  on 
particle 
slip  quantity 
wall  conditions 
freesff'eam  conditions 
initial  conditions 

dimensional 

index  for  dependent  variables:  i  =  1,  2,  3 
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AIR-SAND  BOUNDARY-LAYER  DEVELOPMENT  STARTING  FROM  THE  SHORE  OF  THE 
GULF  OF  AQUABA,  JORDAN,  AS  VIEWED  FROM  EILAT,  ISRAEL  (PHOTOS  BY 
I.  I.  GLASS). 


1.  INTRODUCTION 


Boundary-layer  flows  of  a  dusty  gas  have  been  investigated  using  several 
analytical  methods:  a  series  method  [1-7],  an  integral  method  [8-11],  and  a 
finite-difference  method  [12-15].  All  the  work  mentioned  above,  however, 
dealt  with  incompressible-flow  cases  for  the  gas  phase.  Very  few  authors 
[16-18]  considered  the  problems  of  compressible  boundary-layer  flows  where 
the  density  of  the  gas  phase  can  be  changed  due  to  compressibility.  As 
pointed  out  by  Singleton  [17],  Chiu  [16]  employed  incorrect  boundary-layer 
equations  and  assumed  that  the  particle  density  is  constant.  Singleton 
extended  Marble's  analysis  [1]  to  compressible  boundary-layer  flows.  He 
applied  the  cjordinate  perturbation  method  and  obtained  asymptotic  solutions 
for  two  limiting  regions  (see  Fig.  1):  for  the  large-slip  (or  quasi -frozen) 
region  near  the  leading  edge  (I)  and  for  the  small-slip  (or 
quasi -equi 1 i bri um)  region  far  downstream  of  the  leading  edge  (III).  Zhao 
[18]  used  a  similar  series-expansion  method  and  improved  Singleton's 
analysis.  However,  these  series  solutions  in  the  form  of  asymptotic 
expansions  could  provide  only  one  term  in  addition  to  the  frozen  or 
equilibrium-flow  values,  owing  to  the  complexity  of  the  problem.  Moreover, 
this  solution  does  not  provide  any  information  on  the  boundary-layer 
development  in  the  nonequilibrium  transition  region  where  the  slip  is 
moderate.  A  thorough  understanding  of  the  compressible  dusty-gas 

boundary-layer  over  the  entire-flow  region  is  important,  since  these  flows 
have  practical  applications  in  many  scientific  and  technical  fields  such  as 
solid  rocket  exhaust  nozzles,  nuclear  reactors  with  gas-solid  feeds, 
ablation  cooling,  blast  waves  moving  over  the  Earth's  surface,  conveying  of 
powdered  materials,  fluidized  bed  and  environmental  pollution,  as  mentioned 
in  Refs.  [12,  19]. 

In  the  present  paper,  the  behaviour  of  compressible,  laminar 
boundary-layer  flows  of  a  dusty  gas  over  a  semi -i nf i ni te  flat  plate  along 
the  whole  length  of  the  plate  is  studied  using  a  finite-difference  method. 
The  problem  of  two-phase  suspension  flows  is  solved  in  the  framework  of  a 
model  of  two  interpenetrating  and  interacting  continuous  media,  which  is 
called  a  two-way  coupling  model  or  a  two-fluid  approach  [20,21].  The 
following  assumptions  are  made  in  this  analysis:  (1)  The  gas-particle 

mixture  is  a  dilute  system  where  the  volume  fraction  of  the  particle  phase 
is  neglected.  (2)  The  gas  phase  is  a  perfect  gas.  (3)  The  particles  are 
spheres  of  uniform  size  without  random  kinetic  motion.  There  are  no  mutual 
collisions  or  other  interactions  among  the  particles.  (4)  Only  the  drag  and 
heat-transfer  processes  couple  the  particles  to  the  gas.  The  momentum  and 
energy  exchange  between  the  two  phases  can  be  calculated  from  available 
analytical  solutions  for  the  viscous  flow  field  around  a  single  sphere. 

Finite-difference  methods  of  solution  of  single-phase  boundary-layer 
equations  have  been  studied  for  many  years.  A  review  of  this  work  is  given 
in  Ref.  [22].  Flugge-Lotz  and  Blottner  [23]  developed  an  implicit 
difference  technique.  They  used  a  six-point  scheme  for  the  momentum  and 
energy  equations  and  a  four-point  scheme  for  the  continuity  equation.  This 
finite-difference  procedure  was  applied  successfully  to  various  studies  of 


pure-gas  boundary-layer  flows.  However,  in  the  dusty-gas  case,  the  nature 
of  the  governing  equations  requires  some  changes  which  result  in 
considerable  complexity.  First,  in  addition  to  some  new  interaction  terms 
in  the  conservation  equations  for  the  gas  phase,  there  is  an  extra  set  of 
conservation  equations  for  the  particle  phase.  The  partial  differential 
equations  for  the  gas  phase  are  of  second  order,  while  those  for  the 
particle  phase  are  of  first  order.  Secondly,  there  is  no  correponding  state 
equation  for  the  particle  phase,  since  the  particle  phase  has  no  analog  of 
flow  pressure.  In  order  to  close  the  system  of  basic  equations,  the 
y-momentum  equation  for  the  particle  phase  cannot  be  omitted  as  for  the  gas 
phase.  Finally,  the  flow  properties  of  the  particles  present  quite 
different  features  in  different  flow  regions.  In  the  near  leading-edge 
region,  very  large  velocity  slip  and  temperature  defect  between  the  two 
phases  appear,  whereas  a  quasi-equilibrium  state  can  be  reached  in  the 
far-downstream  region  where  the  flow  profiles  for  the  two  phases  are  almost 
the  same.  Between  these  two  regions,  there  is  a  transition  region  which  is 
characterized  by  a  nonequilibrium  flow.  In  general,  the  two  phases  in  this 
region  have  moderate  differences  in  velocity  and  temperature  across  the 
boundary  layer.  It  is  interesting  to  note  that  in  the  transition  region, 
there  exists  a  special  position  along  the  flat  plate,  which  is  defined  as 
the  critical  point  in  this  analysis.  At  the  critical  point,  the  tangential 
velocity  of  the  particles  at  the  wall  vanishes,  that  is,  there  is  no  slip 
between  the  particles  and  the  gas.  This  is  due  to  the  fact  that  in  the 
two-phase  boundary  layer,  the  gas  decelerates  from  its  freestream  velocity 
at  the  outer  edge  to  zero  at  the  wall  and  then  the  particles  are  retarded  by 
the  gas.  The  velocity  of  the  particles  at  the  wall  may  be  reduced  to  zero 
provided  the  distance  is  long  enough  for  the  particles  to  adjust  to  the  gas. 
Of  course,  equalization  of  the  gas  and  particle  velocities  at  the  wall  does 
not  mean  that  the  disparity  of  the  two  phases  has  died  out  because  across 
the  boundary  layer,  equilibrium  between  the  particles  and  the  gas  is  still 
not  attained.  Nevertheless,  it  can  be  said  that  at  the  critical  point,  the 
two-phase  system  completes  essentially  the  transition  from  a  nonequilibrium 
flow  to  an  equilibrium  flow,  since  the  equilibrium  state  is  reached  first  on 
the  surface  of  the  plate  and  this  process  is  continued  until  the  two  phases 
are  in  equilibrium  across  the  whole  boundary  layer  far  downstream.  As  the 
particles  are  slowed  down,  the  density  of  the  particle  phase  near  the  wall 
increases.  When  the  particle  velocity  becomes  zero,  the  particles  tend  to 
accumulate  at  the  wall.  In  other  words,  deposition  of  the  particles  at  the 
wall  may  occur  if  there  is  no  diffusion.  Therefore,  as  discussed  by  Soo 
[24,25],  there  are  two  possible  situations  when  the  particle  velocity 
decelerates  to  near  zero: 

(1)  For  large  particles,  their  Brownian  motion  is  neglected,  the  particles 
slowed  down  at  the  wall  will  deposit  and  form  a  sliding  layer  (or  bed  of 
particles).  This  compacted  layer  may  build  up  or  erode  away,  and  even  a 
steady  equilibrium  condition  may  be  achieved  when  the  shear  stresses  in 
this  dense  layer  of  the  particulate  matter  and  in  the  suspension  mixture 
are  equalized.  The  velocity  at  which  such  a  layer  moves  depends  on  the 
materials  and  surfaces  of  the  particles  and  the  wall.  Because  of 
deposition  of  the  particle  phase,  the  density  of  the  particles  at  the 


wall  becomes  very  large.  However,  if  the  particle  density  is  too  high, 
the  present  analysis  will  fail  since  the  assumptions  concerning  the 
interactions  between  the  two  phases  or  among  the  particles  in  this  paper 
can  be  considered  correct  just  for  a  dilute  gas-particle  system. 

(2)  For  small  particles,  the  Brownian  diffusion  is  significant  in  the  region 
near  the  wall,  although  the  intensity  of  Brownian  motion  of  the 
particulate  cloud  is  usually  small  across  the  boundary  layer.  The 
density  of  the  particle  phase  at  the  wall  is  then  controlled  by  the 
Brownian  diffusion  process.  The  layer  of  deposited  particles  may  not 
exist  at  all  because  the  diffusion  due  to  the  Brownian  motion  prevents 
the  formation  of  a  dense  bed  of  particles.  It  is  shown  in  Soo's 
analyses  that,  if  the  Schmidt  number  of  the  particle  Brownian  diffusion 
is  of  order  unity  or  less,  the  whole  two-phase  system  behaves  like  a 
gaseous  mixture  and  the  density  profiles  for  the  particle  phase  reduces 
to  its  original  one  as  in  the  freestream.  Soo  studied  only  the 
incompressible  boundary-layer  case.  Marble  [1,26]  treated  the 
case  of  compressible  flows  and  obtained  a  similar  result.  Some  other 
studies  [17,18]  on  the  compressible  boundary-layer  flow  of  a 
gas-particle  mixture  came  to  the  same  conclusion.  From  the  asymptotic 
solution  for  the  small-slip  limit,  it  is  found  that  the  zeroth-order 
approximation  of  the  dimensionless  density  for  the  particle  phase  is  the 
same  as  that  for  the  gas  phase.  This  means  that  the  loading  ratio  of 
the  particles  is  constant  across  the  boundary  layer  and  equal  to  its 
original  value  in  the  freestream.  Physically,  in  this  quasi-equilibrium 
flow  region,  the  particles  always  remain  attached  or  fixed  to  their 
original  gas  mass  and  move  together  with  the  gas.  Then  the  gas-particle 
mixture  behaves  like  a  perfect  gas  with  the  modified  properties.  This 
implies  that  the  flow  process  in  the  small-slip  region  is  mainly 
diffusion-controlled  for  both  the  gas  and  particles.  Therefore,  in  this 
paper,  it  is  assumed  that,  after  the  critical  point,  the  particle 
density  is  determined  from  the  gas  density  and  the  loading  ratio.  Using 
these  considerations,  the  finite-difference  schemes  for  the  dusty-gas 
boundary-layer  flows  can  be  constructed.  In  this  analysis,  the 
f  i  ni  te-di  f  ference  scheme  developed  by  Flugge-Lotz  and  Blottner  [22,23] 
were  employed  for  the  gas  phase.  For  the  particle  phase,  a  four-point 
scheme  was  used.  For  comparison,  the  six-point  scheme  was  used  to  solve 
the  x-momentum  and  energy  equations  of  the  particles,  employing 
additional  boundary  conditions  obtained  from  the  compatibility 
conditions.  After  the  critical  point,  very  simple  compatibility 
conditions  for  the  tangential  velocity  and  temperature  of  the  particles 
can  be  derived:  at  the  wall,  the  particles  have  the  same  velocity  and 
temperature  as  the  gas. 

With  this  finite-difference  scheme,  the  flow  properties  of  the  dusty-gas 
boundary  layer  over  the  entire  length  of  a  semi -infinite  flat-plate  were 
calculated  numerically.  The  flow  profiles  of  u,  v  and  T  for  the  two  phases 
are  presented  at  different  distances  from  the  leading  edge.  From  these 
results,  it  is  shown  that  the  boundary-layer  flows  of  a  dusty  gas  have 
different  characteristics  in  the  three  distinct  regions.  In  the  large-slip 


region,  the  particles  have  a  little  deviation  from  their  freestream  uniform 
motion  and  then  the  differences  in  the  flow  quantities  of  the  two  phases  are 
quite  large.  While  slipping  through  the  gas  downstream,  interaction  between 
the  two  phases  increases  the  gas  velocity  and  temperature  but  decreases  the 
particle  velocity  and  temperature  as  well.  Thus,  in  the  transition  region, 
the  differences  in  the  flow  properties  of  the  two  phases  are  significantly 
reduced.  Of  course,  the  particles  and  the  gas  are  still  in  nonequilibrium. 
In  this  region,  the  velocity  slip  and  temperature  defect  are  moderate 
compared  with  those  in  the  other  two  limiting  regions.  In  the  small-slip 
region  far  downstream,  the  flow  profiles  for  the  particle  phase  become 
almost  identical  with  those  for  the  gas  phase,  that  is,  the  two  phases 

approach  nearly  equilibrium  and  the  slip  quantities  are  very  small.  In 

fact,  the  only  reason  the  particles  do  not  actually  attain  the  local  gas 

velocity  and  temperature  is  that  slip  is  induced  along  the  gas  streamlines 
by  the  gas  retardation  associated  with  thickening  of  the  gas  boundary  layer. 
In  addition  to  the  flow  profiles  of  u,  v  and  T,  some  boundary-layer 

characteristic  quantities  of  interest,  i.e.,  the  shear  stress  and 

heat-transfer  rate  at  the  wall  and  the  displacement  thickness,  are 
calculated  in  this  analysis.  It  is  noted  that  owing  to  the  presence  of  the 
particles,  the  shear  stress  and  heat  transfer  increase  while  the 

displacement  thickness  decreases  in  the  case  of  laminar  boundary-layer 
flows,  since  the  interaction  between  the  gas  and  particles  causes  an 
increase  in  the  gas  velocity  and  temperature. 

In  this  paper,  the  quasi-frozen  flow  properties  in  the  near  leading-edge 
region  and  the  quasi-equilibrium  flow  properties  in  the  far-downstream 
region  were  compared  with  the  corresponding  asymptotic  values  [27],  The 
agreement  was  very  good.  For  the  finite-difference  solution  in  the 

nonequilibrium  transition  region,  it  is  found  that  the  results  are 

physically  reasonable.  Although  it  is  not  possible  at  present  to  make  any 
direct  comparison  between  our  finite-difference  solution  and  other  relevant 
results,  since  there  are  no  experimental  or  other  analytical  data  available 
for  the  nonequilibrium-flow  region.  Nevertheless,  the  fact  that  the 
finite-difference  solution  in  the  far-downstream  region  agrees  quite  well 
with  the  asymptotic  small-slip  solution  provides  confidence  in  the 
difference  solution  for  this  transition  region,  since  the  boundary-layer 
equations  are  parabolic,  which  is  classified  as  a  marching  problem  [30]. 
Thus,  the  solution  procedure  of  finite  difference  begins  with  certain 
initial  profiles  at  or  near  the  leading  edge,  then  through  the  large-slip 
region,  the  transition  region,  and  finally  ends  in  the  small-slip  region 
downstream.  It  is  clear  that  the  finite-difference  solution  for  the 
small-slip  region  would  not  be  correct  if  there  were  some  mistakes  in  the 
difference  solution  for  the  transition  region. 

The  numerical  study  of  boundary-layer  flows  in  dusty  gases  provides  a 
good  introduction  to  the  dynamics  of  a  two-phase  system.  The  quasi-frozen 
flow,  nonequilibrium  flow  and  quasi-equilibrium  flow  are  all  encountered  and 
analysed  using  the  finite-difference  method.  The  difference  solution  gives 
the  complete  and  exact  information  about  modifications  of  the  boundary-layer 


flows  due  to  gas-particle  interaction.  Moreover,  it  provides  a  basis  for 
the  experimental  investigation  of  dusty-gas  boundary-layer  flows. 

2.  MATHEMATICAL  DESCRIPTION  OF  COMPRESSIBLE,  LAMINAR,  DUSTY-GAS  BOUNDARY- 
LAYER  FLOWS 

The  basic  boundary-layer  equations  for  steady,  two-dimensional, 
compressible,  laminar,  dusty-gas  flows  over  a  flat  plate  are  given  by  [27] 

Continuity : 


d  p*g*  +  d  p*v*  =  0 
3x*  by* 


(2.1) 


Momentum: 


P*  (u*  +  v*  M)  =  -A  (p*  &£)  +  pt*(u*  -  u*)  Hr  £  D  (2.2) 

H  ^  a**  av*J  as/*  v  av*;  P  x*  n* 


5y* 


Energy : 


p*cp(u*  f£  *  p>  =  p  <k*  fD  *  “'p)2  *  PP[(UP  '  “*> 


+  (vp*-v*,2]Tf|°*3??ppc5(T5-T*,i^Nu 


U  a, 


(2.3) 


State: 


p*  =  p*R*T* 


(2.4) 


for  the  gas  phase,  and 
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for  the  particle  phase, 
equilibrium  length  \*  is 


The  starred  quantities  in  Eqs.  (2.1)-(2.8)  have  dimensions.  The 
independent  variables  are  the  space  coordinates  x*  and  y*  which  are  parallel 
and  perpendicular  to  the  wall,  respectively.  The  dependent  variables  are 
the  density  p*,  the  velocity  components  u*  and  v*  and  the  temperature  T*  for 
the  gas  phase  as  well  as  the  corresponding  quantities  pi,  u*,  v£  and  T*  for 
the  particle  phase.  For  the  flat-plate  boundary  layer,  the  gas  pressure  p* 
is  constant  and  equal  to  its  freestream  value.  Hence,  in  the  dusty-gas 
boundary-layer  problem,  there  are  eight  simultaneous  equations  with  eight 
unknowns  so  that  this  system  is  closed.  Of  course,  it  is  required  that  the 
other  physical  quantities  appearing  in  Eqs.  (2.1)-(2.8)  are  known  functions 
of  the  flow  variables.  The  normalized  drag  coefficient  D  and  the  Nusselt 
number  Nu  can  be  expressed  in  terms  of  the  slip  Reynolds  number  and  the 
Prandtl  number.  Here,  the  normalized  drag  coefficient  D  is  defined  as 


In  the  equations  above,  the  particle  velocity- 
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where  CD  and  CD  are  the  real  drag  coefficient  for  the  flow  situation  under 
consideration  an9  the  drag  coefficient  from  the  Stokes  relation.  In  this 
analysis,  only  Stokes'  relation  is  used  and  consequently  D  =  1.0  and  Nu  = 
2.0.  Regarding  the  thermodynamic  properties,  the  following  assumptions  are 
made:  (1)  The  specific  heats  for  the  gas  and  particle  phases  (c£  and  c|) 
are  constant;  (2)  the  Prandtl  number  of  the  gas  (Pr)  is  constant;  (3)  the 
viscosity  coefficient  of  the  gas  (n*)  has  a  power-law  form  with  temperature. 
Consequently,  the  expression  for  the  gas  viscosity  is  given  by 


P* 
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where  to  is  the  power  index  for  the  viscosity  coefficient. 

It  is  advantageous  to  write  the  basic  equations  and  relative  expressions 
in  nondimensional  form  before  the  numerical  computations  are  performed.  For 
the  investigation  of  two-phase  boundary-layer  flows,  it  is  convenient  to 
choose  the  velocity-equilibrium  length  \*  as  the  characteristic  length,  and 
the  following  nondimensional  quantities  are  defined. 


where  the  flow  Reynolds  number  based  on  the  particle  velocity-equilibrium 
length  Re*  is 
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Now,  a  nondimensional  form  of  the  boundary-layer  equations  results  in 
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where  the  gas  Eckert  number  Ec,  the  gas  Prandtl  number  Pr  and  the  ratio  of 
specific  heats  for  the  two  phases  a  are  respectively  defined  as 
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The  viscosity  relation  in  nondimensional  form  can  be  written  as 

p  = 


(2.20) 


In  order  to  obtain  a  unique  solution  to  the  problem,  it  is  necessary  to 
satisfy  the  boundary  conditions.  Inspecting  the  basic  equations 
(2.12)-(2.19) ,  there  are  seven  partial  differential  equations  and  two  of 
them  are  of  second  order.  Therefore,  nine  boundary  conditions  should  be 
specified.  If  the  partcle  phase  is  in  equilibrium  with  the  gas  phase  in  the 
external  flow,  the  nondimensional  boundary  conditions  are  given  by: 


(1)  At  the  wall  of  the  flat  plate 


u(x,  0)  =  0, 


v( x ,  0)  =  0,  T(x,  0)  =  Tw, 
Vp(x,0)  =  0 


(2.21) 


(2)  At  the  outer  edge  of  the  boundary  layer 


u(x,  ®)  =  1, 


T( x,  «)  =  1, 


up(x,  »)  =  1, 


Tp(x,  ®)  =  1,  Pp ( x ,  <=)  = 


(2.22) 


where  the  mass  loading  ratio  of  the  particles  p  is 


Besides,  owing  to  the  parabolic  character  of  boundary-layer  equations, 
the  initial  profiles  of  the  dependent  variables  are  required  across  the 
boundary  layer  at  some  point  x0: 


u(x0,  y)  =  u 0{y ) , 

T(x0,  y)  =  T0{y) , 

up(x0,  y)  =  upo(y), 

Tp(x0.  y)  =  TPo(y). 

At  the  initial  position  x0,  the  fini 
and  then  proceeds  downstream. 


v(x0,  y)  =  v 0(y) 

p(x0.  y)  =  p0(y) 

(2.23) 

Vp(x0,  y)  =  vpo(y) 

pp(xo*  y)  =  pPo(y) 

;e-di fference  solution  procedure  starts 


3.  FINITE-DIFFERENCE  SCHEMES  AND  RESULTING  FINITE-DIFFERENCE  EQUATIONS 

The  basic  boundary-layer  equations  (2.12)-(2.19)  with  the  boundary 
conditions  (2.21)-(2.22)  and  the  initial  conditions  (2.23)  can  be  solved 
numerically  using  a  finite-difference  method.  In  this  way,  the 
partial  differential  equations  are  approximated  by  finite-difference 
equations  and  the  flow  field  is  divided  into  a  rectangular  grid  or  mesh. 
Generally,  either  equal  or  unequal  intervals  can  be  used.  In  this  report, 
equal  intervals  in  the  x-direction  and  unequal  intervals  in  the  y-direction 
were  used  in  order  to  reduce  the  computation  time  (see  Fig.  2).  The  step 
size  in  the  y-direction  was  increased  in  a  geometric  progression  as 

*yn  .  , 

A/n-1 

where  K  is  a  constant  and  it  is  set  with  a  value  slightly  greater  than 
unity.  When  K  =  1.0,  the  unequal-interval  mesh  reduces  to  an  equal-interval 
mesh.  In  the  difference  procedure,  it  is  assumed  that  the  flow  quantities 
are  known  at  the  grid  points  in  the  column  (m)  and  unknown  at  the  grid 
points  in  the  column  (m+1).  The  computation  starts  stepwise  downstream  with 
the  initial  profiles. 

When  the  finite-difference  scheme  is  employed,  the  derivatives  are 
replaced  by  difference  quotients.  There  are  numerous  ways  of  constructing 
difference  quotients.  For  the  sake  of  stability,  implicit  schemes,  which 
can  be  six-point  or  four-point,  are  used  in  this  analysis  [28]. 

For  the  momentum  and  energy  equations  of  the  gas  phase,  a  six-point 
difference  scheme  was  used.  With  this  scheme,  six  grid  points  (m,  n-1), 
(m,  n),  (m,  n+1),  (m+1,  n-1),  (m+1,  n)  and  (m+1,  n+1)  are  involved.  Any 
function  w(x,  y)  and  its  derivatives  are  evaluated  at  a  mid-point  (m+0,  n): 
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where  the  weighting  factor  9  can  be  chosen  as  any  value  between  zero  and 
unity.  When  9  =  0.5,  it  reduces  to  the  six-point  Crank-Nicolson  scheme 

where  the  truncation  error  is  of  order  (Ax2)  [29],  When  9  takes  the  value 
of  zero  or  unity,  it  gives  the  full  explicit  or  implicit  scheme, 

respectively.  The  last  two  schemes  involve  only  four  grid  points  and  have  a 
truncation  error  of  order  (Ax).  But  with  respect  to  the  variable  y,  all 
three  schemes  above  are  of  second  order  (Ay2),  since  the  central  difference 
formula  is  used  for  derivatives  in  the  y-direction. 

A  four-point  difference  scheme  is  applied  to  the  gas  continuity 
equation.  In  this  scheme,  four  grid  points  (m,  n-1),  (m,  n),  (m+1,  n-1)  and 

(m+1,  n)  are  included  and  all  the  values  of  the  function  w(x,  y)  and  its 

derivatives  are  calculated  at  a  mid-point  (m+9,  n  -  1/2): 

W  = -1  [0(Wm+i>n  +  +  (l-0)(Wm>n  +  wm,n-l)]  (3.5) 
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When  9  =  0  or  9  =  1  in  the  above  formulae,  the  truncation  error  in  the 
x-direction  is  of  order  (Ax).  When  9  =  0.5,  the  scheme  is  known  as  the 
four-point  Wendroff  scheme  and  the  truncation  error  in  the  x-direction  is  of 
second  order  (Ax2)  [29].  However,  from  experience  in  the  present  analysis, 
the  Wendroff  scheme  may  produce  an  oscillation  in  the  normal  velocity  of  the 


gas  phase.  It  was  found  that  this  oscillation  problem  can  be  avoided  by 
using  0  =  0.75,  which  produces  a  discretization  error  of  order  (ax1*1’).  In 
the  y-direction,  this  four-point  scheme  has  a  truncation  error  of  order 
(Ay2)  as  in  the  six-point  scheme. 

For  the  particle  phase,  due  to  the  stability  requirement,  the 
y-deri vati ves  are  approximated  by  backward  difference  quotients  instead  of 
the  central  difference  quotients  which  are  used  in  the  above  schemes  for  the 
gas.  Then,  another  four-point  difference  scheme  for  the  particles  is 
constructed  as  follows.  The  function  w(x,  y)  and  its  derivatives  are 
estimated  at  the  point  (m+0,  n).  The  derivatives,  both  in  the  x-  and 
y-di rections,  are  replaced  by  backward  quotients. 

For  the  x-momentum,  energy  and  continuity  equations  of  the  particles, 
the  grid  points  (m,  n) ,  (m,  n+1),  (m+1,  n)  and  (m+1,  n+1)  are  involved  and 
the  difference  scheme  is 
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For  the  y-moinentum  equation,  the  grid  points  (m,  n) ,  (m,  n-1),  (m+1,  n) 
and  (m+1,  n-1)  are  involved  and  the  difference  scheme  takes  another  form  for 
the  y-derivative  instead  of  Eq.  (3.10) 
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The  function  w  and  the  x-derivative  dW/dx  have  the  same  forms  as  Eqs.  (3.8) 
and  (3.9). 


Similarly,  when  9  is  equal  to  zero  or  unity,  which  represents, 
respectively,  the  explicit  or  implicit  scheme,  the  above  four-point  schemes 
(3.8)-(3.11)  reduce  to  three-point  schemes.  These  schemes  have  a  truncation 
error  of  first  order  (Ay).  For  stability  consideration,  as  mentioned 
before,  the  value  of  9  is  chosen  to  be  0.75. 


With  the  above  formulae,  the  finite-difference  equations  for 
compressible,  laminar  boundary-layer  flows  of  a  dusty  gas  over  a 
semi-infinite  flat  plate  are  given  by: 


(1)  The  momentum  equation  of  gas  phase, 

%  um+l,n+l  +  3n  um+l,n  +  ^n  um+l,n-l  ”  8n 


(n  =  2,  3, ...N-1)  (3.12) 
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"-F'rrjr  hj  *v»  ,r'tw  rm 


where 


%  =  an^pv  “  p  TyWe.n  "  cntJiiH-9,n 

Bn  =  (P^m+G.n  +  an^2“^^pv  “  p^yW9,n  +  cn^+^pm+9,n 

+  0Ax(^MD)m+1>n 

=  -an^2(pv  “  p  Vnt+9,n  "  cn^pm+9,n 

=  [(Pu)m+9,n  ”  (I- 0) &•(  P^j ^)m,n  Jum, n  “  bn(Pv  "  p  ^y^m+G.n^m.n 


+  ^nPm+G.n A2um,n  +  ^  ^)up^^m+9,n 


(2)  The  energy  equation  of  gas  phase, 

2  2  2  2 
An  ^m+l,n+l  +  Bn  ^m+ 1 , n  +  ^n  ^in+l,n-l  ”  Bn 


(n  =  2,  3....N-1)  (3.13) 


where 


2  1 

An  =  dn(pv  ‘  TyL+9,n  "  cn^U@,n 


Pr 


Bn  ~  (pu^m+9,n  +  an^K2_1HPv  "  -jjrp  Ty  )m+ 9,n  +  cn(K+1)(pr  \n+9,n 
+  9AX^  Pp^Ul.n 

Cn  =  _anK2(Pv  “  Tf-  Ty)m+9,n  '  cn^p:U9,n 


Dn  [(pu)ni+9,n  "  (l-6)Ax(— L_  Pp ^m,n  ^m,n 

bn(pv  -  ii-  Tyj^Q^p  ^^m+9,n  A2^m,n 

+  Ax  {Ec u(Uy ) 2  +  EcPp[(up  -  u)  2  +  _L-  (vp  -  v)2]mD  +  -L  PpTpMNu}m+9tn 
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2 

For  the  coefficient  Dn,  the  term  (vp  -  v)2/ReQ0  is  a  small  quantity  and  can 
be  neglected.  In  all  the  above  coefficient  expressions,  some  parameters  are 
given  by 

a  =  QA*  b  = 

n  (K+l )  Ayn  *  n  (K+l )  Ayn 


29KAx  d  =  2(1- 9)K Ax 

(K+l ) Ayn  2*  n  (K+l )  Ayp  2 


^m.n 

Kn,n+1  + 

(K2-1>V„  - 

^2^m,n-l 

(3.14) 

A2wm,n 

^m,n+l 

(K+1)wm>n  t 

^m,n-l 

(Wy  ^m.n 

_  Hn,n+1  + 

(Kz-l)Wm>n 

-  K2wm,n-1 

(K+l)Ayn 

where  the  function  W  represents  the  flow  properties  such  as  velocity, 
temperature,  density,  etc. 


(3)  The  state  equation  of  the  gas  phase: 

J_  (n  =  1,  2 . N) 


pm+l,n  J 


(3.15) 


m+1  ,n 


(4)  The  continuity  equation  of  the  gas  phase: 

( pvWl  ,n  =  ^  pv^m+l  ,n-l  ~  ^QV^m,n  '  (pv^m,n-lJ 


-  [(pu)m+l,n  "  ^m.n  +  ^m+l.n-l  “  (pu)m,n-l]  (n  =  2»  3’-”N) 

(3.16) 


(5)  The  x-momentum  equation  of  the  particle  phase: 


where 


a!  =  vr 


n  Ayn  Pm+0,n 


=  u 


Pm+9,n  “  ^  VPm+  0,n  +  0Ax( ^m+l ,n 


C3  =  -[(1-8)  ^  v  lu  +  [u  +  (1-8)  Ax  y 
n  Ayn  Pm+0,n  Pm,n+1  Pm+0,n  Ayn  Pm+0,n 


-  (l-0)Ax(pD)rnjn]upm^  +  Ax(upD)m+e>n 


(6)  The  y-momentum  equation  of  the  particle  phase, 

.4  _  4  '*  , 


where 


.  4  „  4  „  4 

An  vn  +  Bn  vn  =  Cn 

n  Pm+l,n  n  Pm+l,n-l  n 


(n  =  2,  3, . . . ,  N) 


4 

An  =  u 


+  9Ax 


Pm+9,n 

%-l 

0Ax 

A*n-1 

pm+e,n 

u 

-  (I-® 

nn+0,n 

% 

-V.,,  * 


’u  „  (1-  8)  Ax  u  0)  Avf  uD)  lv 

P(n+0,n  Ayn_i  Pm+0,n  m,n  Pm,n 

+  [(-l~9)  A*  v  ]vD  +  Ax  ( v  uD ) m+  Q  n 

Ayn_i  pm+0,nJ  Pm,n-1  m  0,n 


(3.18) 


(7)  The  energy  equation  of  the  particle  phase, 
.5  „  5  „  _ 5 


Af1  Vn,n  +  1  +  Bn  TPm+ 1 , n  =  °n 


(n  =  1,  2,  ....  N-l) 


(3.19) 
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where 


\5  =  ®Ax  y 
n  Ayn  Pm+0,n 


Bn  =  un  -  vn  +  8 Ax  (— ?-  pNu  L.i 

n  Pm+e,n  Ayn  Pm+9,n  l3Pr  ■w+l.n 


c;  =  -rlkiilA  v  It  +  ru  +  d-9)&  v 

n  Ayn  PrtH-0,nJ  Pm,n+1  L  pm+9,n  Ayn  Pm+8,n 


-  (l-elAxt^  1*1  U  Vn  +  *(3F  *>•»." 


(8)  The  continuity  equation  of  the  particle  phase, 

6 

j  +  r  —  r 

Pm+l ,n+l 


^n  pp 


+  4  A 


Pm+l ,n  =  C" 


(n  =  1,  2,  ....  N-l) 


where 

a!  =  O'* 


n  Ayn  Pm+0,n 


(3.20) 


B"  *  2*P*1,«  +  ('-29)UPm,n  +  (Vl,»„  -  *W 


+  9U-9)Ax  (v  -  2v  ) 

Ayn  Pm,n+1  Pm,n 


Cn  *  -  »,  *  [(29-l>"0„,  .  * 


Ayn  Pm+9,n  Pm,n+1 


pm+l,n 


m,n 


9(1-Q)  Ax  (y 


Ayn  Pm+l,n+l  Pm+l,n 


-2%..,  J-llziliP*  '  2v_  )]p 


Ay 


Pm,n+1  Pm,n  J  Pm,n 


Using  the  finite-difference  equations  in  the  form  (3.12)-(3.20) ,  a 
stable  and  convergent  numerical  solution  to  the  dusty-gas  boundary-layer 
equations  was  obtained  when  x  is  smaller  than  xcr^ .  After  the  critical 
point  (x  >  *crj),  quite  simple  compatibility  conditions  were  derived  for  the 
particle  velocity  and  temperature.  These  conditions  provided  supplemental 
boundary  conditions  at  the  wall  so  that  the  six-point  scheme  could  be  used 
for  the  x-momentum  and  energy  equations  of  the  particles  when  x  >  xcri.  At 
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the  wall  (y  =0),  with  the  boundary  conditions  u  =  0  and  Vp  =  0,  Eqs.  (2.17) 
and  (2.19)  become 


-(MD)W 


(3.21) 


a 

3Pr 


Tw)UNu)w 


(3.22) 


These  two  equations,  Eqs.  (3.21)  and  (3.22),  are  termed  as  compatibility 
equations  from  which  compatibility  conditions  can  be  derived.  From  Eq. 

(3.21),  it  is  known  that,  as  x  increases,  the  particle  velocity  at  the  wall 
decreases  until  it  becomes  zero.  The  position  of  the  critical  point  is 
determined  by 


upw^xcri )  =  0  (3.23) 

After  the  critical  point  (x  >  xcr^),  the  two  phases  have  zero  velocity  at 
the  wall  and  then  the  drag  vanishes  (Dw  =  0).  Thus,  for  x  >  xcri ,  Eq. 

(3.21)  leads  to 


uD  «  0  (3.24) 

Kw 

Substituting  Eq.  (3.24)  into  Eq.  (3.22)  yields 

TD  »  Tw  (3.25) 

Pw  w 

Equation  (3.25)  is  valid  for  x  >  xc r -j  ,  too.  Equations  (3.24)  and  (3.25) 
mean  that  after  the  critical  point  the  particles  and  gas  are  in  equilibrium 
at  the  wall.  Now,  concerning  the  tangential  velocity  and  temperature  of  the 
particles,  there  exist  two  boundary  conditions:  one  is  at  the  wall  and  the 
other  is  at  the  outer  edge  of  the  boundary  layer,  as  in  the  case  of  the  gas 
phase.  For  the  normal  velocity,  however,  no  such  simple  compatibility 
conditions,  such  as  Eqs.  (3.24)  and  (3.25),  can  be  derived. 


With  the  six-point  scheme,  the  x-momentum  and  energy  equations  of  the 
particle  phase  are  replaced  by  the  following  difference  equations: 


(1)  Momentum: 


A-  Un 

n  Pm+l,n+l 


7  7 

B_  un  +  C_  un 
n  Pm+l,n  n  Pm+l,n-l 


=  D„ 


(n  =  2,  3, 


N-l ) 


(3.26) 


where 


7 

A„  =  a„v„ 
n  n  Pm+0,n 
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B"  ■  Ve.n  +  a"(K2-1)vP«M  +  “,(,flW.'< 


C.  =  -a_K2v. 
n  n  Pm+e,n 


°"  =  lV.,  ■  (1-e)MljD)»>."]Vn  ' 


(2)  Energy: 


A"  TPm+l,n+l  +  Bn  Vn,n  *  C"  TPm+l,n-l  =  (n  =  2*  3 . N_1)(3.27) 


where 


.  8 


An  =  anVPm+0,n 


B„  =  un  +  a_(K2-l)vn  +  9Ax(— —  liNuLi  n 
n  Pm+0,n  n  Pm+0,n  3Pr  ™+1‘n 

Cn  =  -anK2vn 
n  n  Pm+0,n 

of  =  fu„  -  (1-0)  AxLiL  uNu  L  „]T„  -  b_vn  AT„ 

n  Pm+0,n  jPr  m,n  pm,n  n  pm+0,n  Pm 


Therefore,  when  x  >  xcri  »  t,1e  difference  equations  are  composed  of  Eqs. 
(3.12),  (3.13),  (3.15),  (3.16),  (3.18),  (3.26)  and  (3.27)  with  the 

assumption  that  the  particle  density  is  determined  by  pj,  =  Pp.  The  detailed 
derivations  of  all  the  finite-difference  equations  above  are  given  in 
Appendices  A  and  B. 


It  is  noted  that  the  boundary-layer  equations  (2.12)-(2.19)  are  a 
coupled  nonlinear  partial-differential  system.  To  avoid  the  coupling  and 
nonlinearity,  in  the  process  of  discretizing  every  conservation  equation, 
only  one  corresponding  variable  appears  as  an  unknown  in  the  resulting 
difference  equation  and  the  difference  expressions  for  the  products  of  the 
unknown  variables,  functions  or  derivatives  are  chosen  such  that  the  unknown 
variables  appear  linearly  in  the  products.  This  procedure  leads  to  a  linear 
system  of  algebraic  equations  which  are  not  coupled.  As  pointed  out  by 
Blottner  [22]  the  coupling  between  the  equations  results  in  a  tridiagonal 
matrix  which  is  somewhat  more  complicated  to  solve  than  in  the  uncoupled 
case.  In  addition,  for  linear  algebraic  equations,  there  are  several 
effective  methods  of  solution  available.  For  example,  the  Thomas  algorithm 
is  a  very  powerful  and  convenient  technique  to  solve  the  linear  equations 
with  the  tridiagonal  matrix  of  the  coefficients. 
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4.  METHODS  OF  SOLUTION  OF  THE  FINITE-DIFFERENCE  EQUATIONS 


The  methods  of  solution  depend  upon  the  characters  of  the 
finite-difference  equations.  For  the  six-point  scheme,  the  resulting 
difference  equations  constitute  the  system  of  simultaneous  algebraic 
equations  with  a  tridiagonal  matrix  of  the  coefficients.  Advantage  can  be 
taken  of  this  tridiagonal  form  of  the  coefficient  matrix  to  solve  the 
algebraic  equations  by  use  of  the  Thomas  algorithm  [30],  With  the  Thomas 
algorithm,  the  solution  is  obtained  by 


where 


ui  =  Gn  Fn  tyn+l,n+l 
wm+l,n  "  -i - 


(4.1) 


pi  s 
Ln  Hn » 


Fn  =  Bn  -  Cj 


Ln-1 

Fi-i' 


Gn 


-  On  - 


r1 

°n 


5k 

Fn-1 


The  recurrence  relation,  Eq.  (4.1),  can  be  used  to  solve  difference 
equations  (3.12)  and  (3.13),  or  (3.26)  and  (3.27).  Correspondingly,  W1  =  u, 
W2  =  T,  W7  =  up  and  W8  =  Tp. 

By  using  the  following  procedure,  the  flow  profiles  can  be  determined: 

(1)  W/ith  the  boundary  conditions  at  the  wall,  calculate  quantities  E^,  F’, 
Gp  from  the  wall  towards  the  outer  edge. 

(2)  With  the  boundary  conditions  at  the  outer  edge,  calculate  the  flow 
properties  W^  from  the  outer  edge  towards  the  wall. 

After  the  gas  temperature  Tm+}  n  is  known,  the  gas  density  n  can  be 
calculated  directly  from  the  state’  equation  (3.15).  Then,  starting  at  the 
wall  and  using  the  gas  continuity  equation  (3.16),  the  normal  velocity 
vm+i,n  can  be  obtained. 


When  x  <  xcri  ,  the  differential  equations  for  the  particle  phase  are 
discretized  using  the  four-point  difference  schemes.  The  methods  for 
solving  the  difference  equations  (3.17)  to  (3.20)  are  not  difficult.  After 
starting  either  at  the  wall  (only  for  the  y-momentum  equation)  or  at  the 
outer  edge  of  the  boundary  layer,  the  calculations  proceed  consecutively 
from  one  grid  point  to  another  in  recursion,  until  the  whole  boundary  layer 
has  been  traversed: 


i 

Kn+1  ,n 


C1 

’■'n 


W 


i 

m+1 ,n+l 


B 


i 

n 


(i  -  3,  5,  6) 


(4.2) 
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"m+1 ,n 


a  Cn  -  Bn  KUl.n-1 


(i  =  4) 


(4.3) 


where  W3  =  u„  ,  W4  =  v~,  W5  =  Tp  and  W6  =  Pp.  After  the  critical  point 
(x  >  xcri)t  tne  x-momentum  and  energy  equations  of  the  particle  phase  can  be 
discretized  with  the  six-point  scheme  and  then  solved  by  the  Thomas 
algorithm,  as  for  the  gas  phase. 

Before  solving  the  difference  equations  numerically  using  the  methods 
described  above,  some  considerations  are  required: 

(1)  How  to  evaluate  the  coefficient  matrix  elements; 

(2)  How  to  give  the  boundary  conditions  in  a  form  suitable  for  the  numerical 
computation; 


(3)  How  to  obtain  the  initial  profiles; 

(4)  How  to  determine  the  value  of  xcri ,  the  nondimensi onal  coordinate  for 
the  critical  point. 

First,  the  finite-difference  equations  can  be  solved  provided  that  the 
values  of  all  the  coefficients  ,  B^,  and  are  known.  However,  from 
the  expressions  of  these  coefficients,  it  is  seen  that  they  depend  on 
unknown  values  of  the  variables  at  the  grid  line  (m  +  1),  since  the 
difference  scheme  is  an  implicit  one.  This  difficulty  can  be  surmounted  by 
using  an  iteration  procedure.  Of  course,  the  iteration  technique  increases 
the  computation  time  very  much.  The  other  way  to  overcome  the  difficulty  is 
to  use  a  linearization  approximation:  the  quantities  appearing  in  the 

coefficients  are  evaluated  at  the  previous  grid  line  (m)  if  these  quantities 
are  still  unknown  at  the  grid  line  (m  +  1).  Otherwise,  they  take  on  their 
updated  values.  In  this  analysis,  this  linearization  approach  was  employed, 
since  it  is  easier  to  program  and  requires  less  computer  storage.  Of 
course,  it  is  less  accurate  compared  with  the  iteration  approach.  However, 
satisfactory  accuracy  can  be  achieved  by  reducing  the  step  size. 

Second,  in  order  to  solve  the  system  of  simultaneous  algebraic  equations 
at  every  new  grid  line  (m  +  1),  it  is  necessary  to  have  the  boundary 
conditions  in  a  suitable  form.  In  this  analysis  this  is  straightforward, 
since  there  are  no  derivatives  involved  in  the  boundary  conditions  of  the 
Dirichlet  type.  In  the  finite-difference  scheme,  the  boundary  conditions  at 
the  wall,  Eq.  (2.21),  are  written  as 

um+l,l  ”  vm+l,l  ”  ^m+l,n  '  ^w’ 

(4.4) 

vn  =0 

Pm+1,1 


■a '  %  %  \  •„  •.  _  • 


Similarly,  the  boundary  conditions  at  the  outer  edge,  Eq.  (2.22),  are  given 
in  the  form. 


Jm+ 1 , N 


1.  1 


m+1  ,N 


(4.5) 


UPm+l,N  ~  lf  TPm+l,N  ”  PPm+l,N  "  P 


Here,  concerning  the  outer-edge  boundary  conditions,  another  difficulty 
arises:  How  to  select  the  number  N,  the  maximum  value  for  the  number  of 

grid  points  at  the  column  (m+1).  In  other  words,  in  the  computation 
process,  it  is  required  to  know  how  far  to  calculate  the  flow  variables 
across  the  boundary  layer.  In  order  to  guarantee  that  the  value  of  N  used 
represents  the  freestream  condition,  one  can  specify  a  large  number  N,nax  for 
the  grid  points  at  the  last  grid  line  (mmax)  far  downstream,  where  the 
computation  terminates.  For  all  the  previous  grid  lines  (0  <  rn  <  mmax),  the 
same  number  N(nax  is  used  to  define  the  outer  edge  of  the  boundary  layer. 

This  .method  is  direct  but  inefficient,  since  it  needs  more  computation  time. 

There  is  another  approach  in  which  a  special  value  of  N  is  chosen  for  a 

given  grid  line  (m) .  In  the  latter  method,  N  is  determined  as  follows 

[31]: 


(1)  It  is  assumed  that  at  the  new  grid  line  (m  +  1)  is  equal  to  Nm,  the 

number  of  grid  points  at  the  previous  line  (m). 

(2)  The  finite-difference  equations  are  solved  with  the  assumed  and  the 

values  of  the  flow  quantities  at  the  last  two  consecutive  grid  points, 

Wm+1,N-1  and  Wm+1,N’  can  be  found* 

(3)  The  difference  between  Wm  +  j  and  Wm+i  N  is  compared  with  a  certain 
small  quantity  e.  If  the  condition  of  smooth  conjugation 


lWm+l,N  ‘  Wm+1,N-1 


(4.6) 


is  satisfied,  the  selection  of  is  correct  and  the  computation  can 

proceed  to  the  next  step. 


(4)  If  the  condition  (4.6)  is  not  fulfilled,  it  is  required  to  assume  = 

Nn  +  1  and  then  to  obtain  the  new  values  of  Wm+j  ^  and  Wm+^  If  the 
condition  (4.6)  is  not  fulfilled  again,  it  is  necessary  to  increase  Nm+^ 
with  unity,  and  so  forth,  until  the  smooth-conjugation  condition  is 
satisfied.  With  this  method,  the  number  of  grid  points  across  the 
boundary  layer  varies  as  the  thickness  of  the  boundary  layer  increases. 

Next,  in  order  to  initiate  the  computation,  the  initial  flow  profiles 
must  be  given.  In  most  pure-gas  boundary-layer  studies,  the  initial 
profiles  are  obtained  from  similarity  solutions.  For  dusty-gas 
boundary-1 ayer  flows,  however,  no  analogous  similarity  solutions  exist.  In 
previous  work  on  finite-difference  solutions  for  incompressible  dusty-gas 
boundary-1 ayer  equations,  the  initial  profiles  were  specified  in  two  ways: 
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(1)  The  Blasius  similarity  profiles  were  chosen  for  the  gas  phase  and 
uniform  profiles  for  the  particle  phase  [14]. 

(2)  The  initial  profiles  were  obtained  by  using  an  integral  method  [13, 
15]. 

It  is  well  known  that  all  the  integral  methods  for  boundary-layer  analysis 
do  not  attempt  to  satisfy  the  basic  equations  at  every  point;  instead,  they 
guess  or  assume  a  suitable  expression  for  the  velocity  and  temperature 
profiles  and  satisfy  the  boundary-1 ayer  equations  only  on  an  average 
extended  over  the  thickness  of  the  boundary  layer.  In  general,  the  initial 
profiles  obtained  from  integral  methods  are  quite  approximate.  From  the 
studies  of  the  behaviour  of  dusty-gas  boundary-layer  flows  near  the  leading 
edge  [17,  18,  27],  it  is  also  known  that  the  similarity  profiles  for  the  gas 
phase  and  the  uniform  profiles  for  the  particle  phase  are  the  zeroth-order 
approximation  in  the  large-slip  region.  The  zeroth-order  asymptotic 
profiles  are  physically  reasonable  and  were  tested  in  this  analysis.  More 
accurate  initial  profiles  up  to  the  first  order  can  be  obtained  from  the 
asymptotic  large-slip  solution  to  the  compressible  laminar  boundary-layer 
equations  for  the  gas-particle  flow  over  a  semi-infinite  flat  plate  by  using 
a  series-expansion  method  [27].  Thus,  it  is  suggested  here  to  employ  the 
first-order  asymptotic  solution  as  the  initial  profiles.  However,  in  this 
approach,  it  is  required  to  obtain  the  asymptotic  solution  first  and  then  to 
solve  the  difference  equations  starting  at  a  given  initial  position  which 
may  be  very  near  the  leading  edge  but  cannot  be  exactly  at  the  leading  edge, 
since  the  asymptotic  solution  involves  a  singularity  at  the  leading  edge. 
Wu  [32]  once  proposed  the  following  type  of  initial  profiles  at  the 
leading  edge  in  pure-gas  cases: 

(1)  The  tangential  velocity  u*  and  temperature  T*  have  their  freestream 
values  at  all  the  grid  points  across  the  boundary  layer  except  at  the 
wall. 


(2)  At  the  wall,  the  tangential  velocity  u*  is  zero  and  the  temperature  T* 
corresponds  to  the  wall  temperature. 

(3)  The  normal  velocity  v*  is  assumed  zero  at  all  the  grid  points. 

Clearly,  this  method  is  very  advantageous  for  starting  numerical  computation 
without  any  preliminary  calculations  of  initial  profiles.  Flugge-Lotz  and 
Blottner  [23]  studied  the  possibility  of  using  the  Wu-type  initial  profiles 
in  pure-gas  boundary-1 ayer  cases  and  concluded  that  the  Wu-type  of  initial 
profiles  can  give  reasonable  results  if  the  proper  mesh  sizes  are  chosen. 
In  the  dusty-gas  case,  similar  initial  profiles  can  be  set  up  by  using 

Wu-type  profiles  for  the  gas  phase  and  uniform  profiles  for  the  particle 

phase.  Soo  [24]  applied  such  initial  profiles  to  his  analysis  of 

incompressible,  laminar  boundary-layer  flow  of  a  dilute  particulate 
suspension.  The  profiles  are  termed  here  as  the  extended  Wu-type.  In  this 
report,  the  above  three  different  types,  i.e.,  the  first-order 


asymptotic-type,  the  zeroth-order  asymptotic-type  and  the  extended  Wu-type, 
were  used  respectively  as  the  initial  profiles  in  order  to  compare  them. 


As  mentioned  before,  the  asymptotic  types  of  initial  profiles  can  be 
obtained  from  the  large-slip  solution  [27].  However,  because  of  the 
different  notations,  certain  relations  must  be  established  with  the 
asymptotic  large-slip  solution  as  follows: 


x  =  (x) 


asy  * 


y  =  U2x  ti) 


asy 


u  -  (u)asy »  ^  ~  f  ( TP  -  0  )asv >  ^  “  (T)a<;v»  P  ~  (  p) 


/2x 


asy1 


asy ! 


'asy 


(4.7) 


u  =  (_1^£)  ,  vn  =  -(^x  (ft  +  i  -  _D  ^!E)) 

p  pp  bn  asy  p  pp  ox  2x  2x  bn  ]  asy 


Tp  *  (Tp)  .  Pp  =  P(Pp) 
p  p  asy  p  p  asy 


Similarly,  a  set  of  relations  can  be  written  for  connecting  the 
finite-difference  solution  with  the  asymptotic  small-slip  solution.  In  this 
analysis,  nondimensi onal  slip  quantities  for  the  particle  phase  are  defined 
as 


Us  =  Up  -  u. 


(4.8) 


The  corresponding  relations  are 


x  =  (x) 


asy’ 


y  =  J_  ( /2x  n) 
/1+p 


asy 


us  =  (us) 


asy 


vs  = 


1 


(-L)  .  T,  =  (TJ 


/1+p  /2x  asy 


asy 


(4.9) 


These  expressions  are  useful  when  comparing  the  finite-difference  solution 
with  the  asymptotic  small-slip  solution.  In  Eqs.  (4.7)  and  (4.9),  the 
subscript  'asy'  denotes  the  asymptotic  solution. 


Finally,  as  pointed  out  earlier,  after  the  critical  point  (x  >  xcr-j)>  it 
is  assumed  that  the  particle  density  is  equal  to  the  gas  density  times  the 
mass  loading  ratio  of  the  particles.  It  is  equivalent  to  assume  that  there 
is  no  accumulation  of  particles  on  the  surface  of  the  plate  and  the  flow  is 
then  mainly  diffusion-controlled  for  the  particle  phase  as  well  as  for  the 


gas  phase.  In  addition,  after  the  critical  point,  quite  simple 
compatibility  conditions  (3.24)  and  (3.25)  are  valid.  The  value  of  xcr1  can 
be  determined  from  the  compatibility  equation  (3.21)  and  the  condition 
(3.23).  Equation  (3.21)  is  an  ordinary  differential  equation  and  the 
solution  upwU)  can  be  obtained  numerically  or  analytically.  For  instance, 

in  the  case  of  the  Stokes  relation  (D  =  1.0),  Eq.  (3.21)  can  be  integrated 
analytically  as 

upw(x)  =  1  -  Mv,X  (4.10) 

where  un  (0)  =  1  is  taken  as  the  initial  condition  at  x  =  0.  From  the 
condition  Up  (xcri-)  =  0,  Eq.  (3.23),  the  critical  value  xcri-  can  be 
determined,  say,  for  the  Stokes  case 

*cri“—  (4-n> 

when  Tw  =  0.5  and  w  =  0.5,  xcr^  =  /2  or  x*ri  =  /2  \*.  If  p*  =  2.5  g/cm3, 
d*  =  10  m,  T*  =  300  K  (or  p*  =  1. 80*10-^  NS/m2)  and  u*  =  500  m/S,  the 
typical  values  of  relaxation  parameters  are  obtained  as:  X*  =  0.386m  and 
x*n-  =  0.546  m. 


5.  RELATIONS  FOR  SHEAR  STRESS,  HEAT  TRANSFER,  AND  DISPLACEMENT  THICKNESS 

Once  the  gas  flow  profiles  across  the  boundary  layer  are  determined, 
some  boundary-layer  characteristics  of  practical  interest  can  be  calculated. 
The  important  quantities  describing  the  behaviour  of  boundary-layer  flows 
are  shear  stress  at  the  wall,  heat-transfer  rate  at  the  wall  and 
displacement  thickness.  They  are  given  in  dimensional  form  as 

(1)  Shear  stress  at  the  wall: 

’w  *  14©  (5.1) 

"  "  dy*  w 


(2)  Heat-transfer  rate  at  the  wall: 


6*  =  /  (1 
o 


(5.2) 


(5.3) 
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The  corresponding  nondiinensi  onal  characteristic  quantities  are  defined  as 


— dte.  (5.4) 

Pt  u*  2 


qw  =  --w--  /Re,  (5.5) 

P®  u*  3 


5  =  £  /Re.  (5.6) 

'"CO 


Substituting  the  nondimensi onal  transformations  (2.11)  into  the  expressions 
(5.1)-(5.3),  the  nondimensional  boundary-layer  characteristics,  Eqs. 
(5.4)-(5.6),  can  be  written  as 


s,  =  *,(£) 

"  by  w 

(5.7) 

i.  .  -  ‘V  (*) 

w  Pr  Ec  by  w 

(5.8) 

ao 

6  =  /  (1  -  pu)dy 

0 

(5.9) 

For  numerical  computations,  it  is  necessary 
relations  in  finite-difference  form.  By  means  of 
gas  velocity  u  and  temperature  T  near  the  wall 
sifficient  accuracy  as 

to  express  the  above 
polynomial  fitting,  the 
may  be  expressed  with 

u  =  au  +  V  +  V2  +  V3 

(5.10) 

T  =  aT  +  byy  +  Cyy 2  +  dyy3 

(5.11) 

Taking  the  derivatives  of  the  above  variables  with  respect  to  y  and  setting 
y  =  0,  the  formulae  for  shear  stress  and  heat-transfer  rate  at  the  wall  are 
obtained  as 


\  "  Mwbu*  %  pr\  bT 

The  values  of  b  and  by  can  be  determined  by  evaluating  Eqs.  (5.10)  and 
(5.11)  at  the  four  grid  points  nearest  the  wall  and  solving  the  resulting 
equations.  Then  the  shear  stress  and  heat-transfer  rate  can  be  given  by  the 
following  expressions  (see  Appendix  C): 


\  = 


(5.12) 


(K2  K  +  1)  [  _  U3  +  U4 

K2Ay1  K^K  +  V  K(K 2  +  K  +  1) 2 


JV_  (K2  +  K  *  1)  [(T  .  x  )  .  T3  -_T1 

>r  Ec  K2iyi  12  1  K(K  *  1 


qw  =  - 


t4  -  Ti  ] 


V  K(K  2  +  K  +  1)  2 


(5.13) 


where  the  subscripts  1,  2,  3,  4  denote  the  four  grid  points  nearest  the  wall 
and  uL  =  uw  =  0,  Tj  =  Tw. 


To  calculate  the  nondiinensional  displacement  thickness  6,  a  three-point 
difference  formula  of  integration  was  used.  The  formula  can  be  applied  to  a 
nonequidi stant  step  size  [33] 


%-l 


r3K  +  2  p 


+  3K  +  1 


K(K  +  1) 


(5.14) 


where  F  =  1  -  pu,  and  the  subscripts  n-1,  n  and  n-t-1  represent  three 
consecutive  grid  points  at  a  given  section. 


6.  COMPUTER  PROGRAM 

The  computer  program  FD8LEP  for  solving  boundary-layer  equations  of  a 
dusty  gas  over  a  semi -inf inite  flat  plate  was  written  in  Fortran  language  on 
a  Perkin-Elmer  computer  system  at  IJTIAS  (see  Appendix  0). 

Only  the  main  features  of  the  calculation  procedure  will  be  reviewed 
here.  A  rectangular-grid  system  indicated  in  Fig.  2  was  adopted,  with  the 
m-lines  in  the  y-direction  (i.e.,  normal  to  the  plate)  and  the  n-lines  in 
the  x-direction  (i.e.,  parallel  to  the  plate).  The  y-axis  is  designated  as 
the  initial  line  for  the  m-lines  and  the  x-axis  for  the  n-lines.  In  the 
finite-difference  procedure,  the  flow  profiles  along  some  m-line  (say,  the 
initial  line  m  =  0)  are  known  and  the  flow  parameters  along  the  (m+l)-line 
have  to  be  determined. 

Given  below  are  the  main  steps  in  the  computation  procedure: 

(1)  Compute  i^+1>n. 

Using  the  known  solution  at  the  m-line  (linearized  conditions)  and  the 
boundary  conditions  on  the  (m+l)-line,  the  new  tangential  velocity  of 
the  gas  at  all  the  grid  points  of  the  (m+1)  line,  um+i  n,  are 
calculated. 

(2)  Test  for  the  outer  edge  of  the  boundary  layer  while  computing  u.  After 

the  boundary  layer  has  been  traversed,  the  two  consecutive  values  of  u, 
u  +i  and  are  compared  to  see  if  the  difference  between  them 

is  fess  than  some  small  positive  number  e.  The  value  of  c  is 
determined  by  the  desired  accuracy  of  the  computation. 
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(3)  Compute  Tm+1>rr 


Using  the  new  values  of  uin+j  with  the  other  linearized  conditions, 
the  new  gas  temperatures  T^+1  n  across  the  boundary  layer  are 
calculated. 


(4)  Compute  Pm+1>n. 

Using  the  new  values  of  Tm+i  n,  the  new  gas  density  profile  is 
calculated  directly. 


(5)  Compute  vm+1>n. 

Using  the  new  values  of  uin+j  n  and  pm+]_  n,  the  new  normal  velocity  of 
the  gas  vm+^  n  across  the  boundary  layer  are  calculated. 

(6)  Compute  u^ 

Pm+l,n 

Using  the  new  tangential  velocity  profile  um+j^  n  for  the  gas  phase  and 

the  linearized  velocity  profiles  u„  and  vn  for  the  particle 

pm,n  pm,n 

phase,  the  new  tangential  velocity  profile  of  the  particles  is 
computed . 


(7)  Compute  vn 

pm+l ,n 

Using  the  new  values  of  vm+3>n  and  up  ^  ^  with  the  linearized  values 

of  v  ,  the  new  normal  velocity  profile  for  the  particle  phase  is 
Hm  ,n 
computed. 


(8)  Compute  T_ 

Pm+1 ,n 

Using  the  new  values  of  u„  ,  vn  and  Tm+i  n  with  the  values  of 

pm+l,n  pm+l ,n  m+1’n 

Tpat  the  previous  line  (m) ,  the  new  temperature  profile  for  the 
particle  phase  is  computed. 


(9)  Compute  pn 

1  m+1  ,n 

Using  the  new  values  of  un  and  v_  with  the  values  of  un,  v_, 

y  Pm+l,n  pm+l,n  p  p 

Pp  at  the  previous  line  (m)  ,  the  new  particle  density  profile  is 

computed.  However,  when  x  >  *Cri»  new  Part’c^e  density  profile 
can  be  obtained  by  setting  pp  i>n. 

(10)  Compute  xw. 


Using  the  first  four  values  of  velocity  for  the  gas  phase  nearest  the 
wall,  i.e.,  um+1  1?  um+1>2,  and  um+1>4,  the  shear  stress  at  the 
wall  is  calculated. 
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(11)  Compute  q  . 

Using  the  values  of  gas  temperature  at  the  four  grid  points  nearest  the 
wall,  i.e.,  Tin+ljl,  T^+i^.  Tm+1,3  and  Tm+1  4>  the  heat-transfer  rate 
at  the  wall  is  calculated.’ 

(12)  Compute  6. 

Using  the  density  and  tangential  velocity  profiles  for  the  gas  phase, 
pm+l  n  and  um+l  n*  ttle  displacement  thickness  can  be  obtained  by 
integration. 

To  advance  the  computation  from  the  (m+l)-line  to  the  (m+2)-line  and  so  on, 
the  same  procedure  (1)-(12)  is  repeated  until  the  desired  value  of  x  is 
reached.  The  flow  diagram  for  the  basic  computer  program  is  shown  in  Fig. 
3.  It  is  noted  that  the  order  of  solving  the  difference  equations  is 
important.  The  equations  for  the  gas  phase  are  solved  first  where  the 
momentum  and  energy  equations  must  be  solved  before  the  continuity  equation. 
Then  the  equations  for  the  particle  phase  are  solved  and  the  energy  and 
continuity  equations  must  be  solved  after  the  momentum  equations. 


7.  NUMERICAL  RESULTS  AND  DISCUSSIONS 


The  present  finite-difference  technique  was  used  to  solve  the 
compressible,  laminar  boundary-layer  flows  of  a  dusty  gas  over  a 
semi-infinite  flat  plate.  The  difference  solutions  for  the  three  flow 

regions  were  obtained:  the  quasi-frozen  flow  region  near  the  leading  edge, 
the  nonequilibrium  flow  region  and  the  quasi -equi 1 i bri urn  flow  region  far 

downstream.  The  asymptotic  solutions  in  the  two  limiting  regions  (the 
large-slip  and  small-slip  regions,  respectively)  were  also  solved 
numerically  in  order  to  independently  verify  the  validity  of  the  present 
implicit  finite-difference  scheme  when  it  is  applied  to  a  gas-particle 
mi xture. 

In  this  analysis,  the  following  set  of  parameters  was  chosen  so  that  the 
finite-difference  results  could  be  compared  with  the  asymptotic  solutions  by 
Singleton  [17]: 

(1)  The  power  index  of  the  gas  viscosity  is  equal  to  0.5  ( to  =  0.5). 

(2)  The  Prandtl  number  of  the  gas  is  equal  to  unity  (Pr  =  1.0). 

(3)  The  Eckert  number  of  the  gas  is  equal  to  unity  (Ec  =  1.0). 

(4)  The  ratio  of  specific  heats  for  the  two  phases  is  equal  to  unity  (a  = 

1.0). 
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(5)  [he  mass  loading  ratio  of  the  particles  is  equal  to  unity  ( p  =  1.0). 

(6)  Stokes'  relation  applies  for  the  interaction  between  the  two  phases  (D  = 

1.0  and  Nu  =  2.0). 

(/)  The  nondi mens i onal  temperature  at  the  wall  is  equal  to  0.5  (T  =  0.5). 

The  flow  in  the  large-slip  region  was  considered  first.  The  initial 
profiles  at  xQ  =  0.005  were  obtained  from  the  first-order  asymptotic 
solution.  The  computation  proceeded  from  xQ  =  0.005  to  x  =  0.505  with  inesh 
parameters  Ax  -  0.001,  Ay  1  =  0.01  and  K  =  1.05.  The  flow  profiles  for  the 

two  phases  at  x  =  0.055,  0.105  and  0.305  are  plotted  in  Figs.  4  to  6.  As 

these  results  show,  there  is  a  very  large  slip  between  the  particles  and  the 
gas  in  this  near-1 eading-edge  region.  Then  the  flow  is  quasi -frozen.  This 
situation  can  be  explained  as  follows.  In  the  freestream,  the  gas  and 
particles  are  in  equilibrium,  that  is,  they  have  the  same  tangential 
velocity  and  temperature  while  their  normal  velocity  is  equal  to  zero.  At 
the  leading  edge,  due  to  viscous  effects,  the  tangential  gas  velocity 
decreases  from  its  freestream  value  at  the  outer  edge  to  zero  at  the  wall 
and  the  gas  temperature  also  changes  from  its  freestream  value  at  the  outer 
edge  to  the  wall  temperature  at  the  wall,  whereas  the  normal  gas  velocity 
acquires  quite  a  large  value.  The  particles,  however,  cannot  accommodate 

these  rapid  changes  and  tend  to  keep  their  original  state  of  motion  in  the 
freestream.  It  takes  some  time  for  the  particles  to  attain  their 

equilibrium  with  the  gas.  The  relaxation  process  of  the  particles  occurs 
throughout  the  ve 1 oc i ty-equ i 1 i bri urn  length  This  two-phase  slip 

phenomenon  implies  that  the  viscous  relaxation  length  for  the  gas  is  much 
shorter  than  the  relaxation  length  for  the  particles  owing  to  the  drag  and 
heat-transfer  interactions  between  the  two  phases.  By  comparing  the 
finite-difference  solution  with  the  asymptotic  solution  in  the  large-slip 
limit,  it  is  seen  that  excellent  agreement  is  obtained  when  x  <  0.1. 

Therefore,  the  asymptotic  large-slip  solution  is  valid  when  x  <  0.1. 

Next,  the  dusty-gas  boundary-1 ayer  flow  in  the  nonequilibrium  region  was 
studied.  The  first-order  asymptotic  profiles  at  xQ  =  0.05  were  taken  as  the 
initial  profiles.  The  grid  parameters  were  Ax  =  0.001,  Ay^  =  0.03  and  K  = 
1.05.  The  flow  profiles  of  the  two  phases  at  x  =  0.55,  1.05,  2.05  and  5.05 

are  shown  in  Figs.  7  to  10,  respectively.  The  numerical  results  indicate 

that  the  slip  between  the  particles  and  the  gas  diminishes  gradually  as  x 
increases.  However,  in  this  transition  region,  the  particles  still  have 
moderate  slip  against  the  gas  and  then  the  two-phase  flow  is  characteri  zed 
by  nonequilibrium.  At  the  critical  point  (x  =  xcr-j),  the  particles  achieve 
equilibrium  with  the  gas  at  the  wall.  After  that  point,  the  two  phases  are 
in  a  state  of  near-equilibrium.  From  the  experience  of  this  analysis,  after 
the  critical  point,  the  four-point  scheme  for  the  particle  continuity 
equation  became  unstable.  It  was  necessary  to  seek  an  appropriate  treatment 
of  the  particle  density.  As  mentioned  before,  based  on  the  fact  that  the 
gas-particle  mixture  acts  like  a  perfect  gas  with  the  total  density  (l  +  j3)p* 
in  the  small-slip  limit,  it  is  assumed  that  the  particle  phase  has  the  local 
density  pjj  =  Bp*,  which  means  that  the  mass  loading  ratio  of  the  particles 
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is  constant  across  the  boundary  layer.  Of  course,  this  assumption 

concerning  the  particle  density  for  x  >  xcn-  may  cause  some  inaccuracy  in 
the  prediction  of  the  flow  properties  in  the  nonequilibrium  region  since 
this  density  distribution  represents  the  zeroth-order  approximation  in  the 
small-slip  limit.  Fortunately,  a  near-equilibrium  state  between  the  two 
phases  is  essentially  reached  after  the  critical  point,  as  the  difference 
solutions  show.  The  effect  of  the  particle  density  on  the  flow  properties 
takes  place  only  through  the  interaction  terms.  Under  the  near-equilibrium 
condition,  these  interaction  terms  will  become  small  of  second  order 
compared  with  the  viscous  terms.  Therefore,  this  approximate  treatment  of 
the  particle  density  is  acceptable.  It  provides  an  approach  to  solve  the 
dusty-gas  boundary-layer  flow  in  the  region  x  >  xcr^.  More  discussion  about 
the  particle  density  is  given  in  Appendix  E.  In  Figs.  7  to  10,  it  is  also 
noted  that  the  relaxation  of  the  tangential  velocities  of  the  two  phases  is 
terminated  effectively  at  about  x  =  5.0  as  well  as  the  temperatures.  In 
contrast  with  tangential  velocity,  there  is  still  an  apparent  difference 

between  the  normal  velocities  of  the  two  phases  up  to  x  ~  20.0.  It  seems  to 
mean  that  the  relaxation  of  Vp  and  v  occurs  over  a  greater  length  than  that 

for  Up  and  u.  In  fact,  the  main  reason  is  that  only  one  mechanism,  i.e., 

the  interaction  between  the  phases,  acts  in  the  relaxation  process  for  up 
and  u  while  two  mechanisms,  the  interaction  between  the  phases  and  the 
continuity  requirement,  both  play  an  important  role  in  the  process  for  vp 
and  v. 

For  the  finite-difference  solution  in  the  small-slip  region,  the 
computation  was  started  at  xQ  =  5.05  and  continued  until  x  =  20.05.  The 
initial  profiles  were  obtained  from  the  finite-difference  solution  in  the 
nonequilibrium  region.  All  the  mesh  parameters  used  in  this  calculation 
were  the  same  as  those  in  the  transition  region.  The  numerical  results  for 
x  =  10.05,  15.05  and  20.05  are  shown  in  Figs.  11  to  13.  It  is  seen  that  the 
quasi -equi 1 i brium  state  between  the  particles  and  the  gas  has  already  been 
reached.  In  Fig.  14,  the  particle  slip  quantities  us ,  v$  and  Ts  obtained 
from  the  difference  solution  are  compared  with  those  from  the  asymptotic 
small-slip  solution.  There  is  very  good  agreement  between  these  two 
solutions.  The  comparison  between  the  finite-difference  solution  and  the 
asymptotic  series-expansion  solution  in  the  small-slip  region  indicates  that 
the  finite-difference  scheme  presented  in  this  paper  has  proved  to  be  a 
useful  method  for  studying  dusty-gas  boundary-1 ayer  flows  and  that  the 
asymptotic  small-slip  solution  is  valid  when  x  >  10. 

It  is  found  from  Figs.  4  to  13  that  the  tangential  velocity  of  the 
particles  in  the  boundary  layer  is  always  greater  than  that  of  the  gas  and 
it  decreases  monotoni cal ly  from  its  freestream  value  to  the  value  at 
equilibrium  with  the  gas  as  one  advances  downstream  from  the  leading  edge. 
Regarding  the  normal  velocity,  however,  the  situation  is  a  little  different. 
Near  the  leading  edge,  the  normal  velocity  of  the  particles  is  smaller  than 
that  of  the  gas,  i.e.,  vp  <  v.  With  increasing  distance  x,  the  normal 
velocity  of  the  particles  becomes  greater  than  that  of  the  gas,  which 
happens  first  near  the  wall  and  then  extends  over  the  whole  thickness  of  the 
boundary  layer.  However,  with  increasing  x,  say  x  =  20,  the  normal 
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velocities  of  the  particles  and  the  gas  reach  an  equilibrium  value  as  for 
the  tangential  velocities.  This  difference  can  be  explained  as  follows. 
The  normal  velocities  for  the  two  phases  are  equal  to  zero  in  the  freestream 
and  induced  to  some  values  in  the  boundary  layer  whereas  the  tangential 
velocities  for  the  two  phases  are  equal  to  the  freestream  value  and  decrease 
in  the  boundary  layer.  In  the  boundary  layer,  the  gas  tangential  velocity 

at  the  wall  vanishes  and  its  distribution  across  the  boundary  layer  is 

similar  to  the  profile  for  the  pure-gas  case.  At  the  leading  edge,  the 

particles  tend  to  keep  their  motion  in  the  freestream  and  consequently  the 
particle  tangential  velocity  is  greater  than  the  gas  velocity.  Then,  owing 
to  the  slip  between  the  two  phases,  a  drag  force  arises  and  the  particles 
decrease  their  tangential  velocity  while  the  gas  increases  its  tangential 
velocity.  Marching  downstream,  the  difference  between  the  tangential 
velocities  for  the  two  phases  becomes  smaller  and  smaller.  With  the 

particle  slip  velocity  approaching  zero,  the  drag  exerted  on  the  particles 
by  the  gas  approaches  zero  as  well.  When  the  particle  velocity  becomes  the 
same  as  that  of  the  gas,  the  interaction  between  the  particles  and  the  gas 
disappears.  Therefore,  the  particle  tangential  velocity  may  become  nearly 
the  same  as  the  gas  tangential  velocity  but  it  cannot  be  smaller  than  that. 
By  contrast,  after  entering  the  boundary  layer  at  the  leading  edge,  the 

normal  gas  velocity  is  induced  first  to  satisfy  the  continuity  requirement 
and  the  normal  velocity  of  the  particles  is  then  induced  due  to  the  drag 
force  exerted  by  the  gas.  Of  course,  this  induced  velocity  for  the 

particles  cannot  be  greater  than  that  of  the  gas.  However,  in  addition  to 
the  interaction  between  the  two  phases,  the  continuity  requirement  for  the 
particle  phase  is  another  important  factor  which  determines  the  changes  in 
the  normal  velocity.  Since  the  tangential  particle  velocity  decreases  in 
the  x-direction,  a  normal  velocity  must  be  induced  to  ensure  mass 
conservation.  Thus,  the  normal  particle  velocity  may  exceed  the  normal  gas 
velocity.  Especially  in  the  region  near  the  wall,  where  the  particle  slip 
velocity  is  quite  large,  the  retardation  of  the  particles  is  considerable 
and  it  results  in  a  rapid  increase  of  the  normal  particle  velocity.  In  the 

region  far  downstream,  where  the  thickness  of  the  boundary  layer  varies  very 

slowly,  the  effect  of  the  interaction  between  the  two  phases  becomes 
predominant  and  the  normal  velocities  of  the  two  phases  tend  to  approach 
each  other  like  the  tangential  velocities. 

The  numerical  results  shown  above  were  obtained  with  the  six-point 

scheme  for  the  x-momentum  and  energy  equations  for  the  particle  phase  after 
the  critical  point.  Instead,  if  the  four-point  scheme  is  still  employed  for 
the  particle  x-momentum  and  energy  equations  after  the  critical  point, 
regardless  of  the  compatibility  conditions,  the  corresponding  results  are 
shown  in  Figs.  15  to  17.  It  is  found  that  some  small  oscillations  appear  in 
the  particle  temperature  profiles  in  the  near-equi 1 ibrium  flow  region,  for 
example,  at  x  =  4.05  and  6.05  (see  Figs.  15  and  16).  The  reason  is  mainly 
that  the  four-point  scheme  has  a  first-order  truncation  error  in  the 

y-direction,  0(Ay),  while  the  six-point  scheme  is  of  second  order,  0(Ay2). 
In  other  words,  the  latter  is  more  accurate  and  it  can  lead  to  a  better 
result.  It  is  also  interesting  to  note  that  these  oscillations  are  bound 

and  damp  out  as  x  increases.  In  Fig.  17  (x  =  10.05),  it  is  seen  that  these 


oscillations  disappear.  It  means  that  in  the  small-slip  region  (x  >  10), 
both  the  six-point  scheme  and  the  four-point  scheme  can  be  used. 

The  finite-di fference  computations  were  made  with  the  mass  loading  ratio 
P  =  0,  using  the  same  difference  scheme  for  the  dusty-gas  boundary  layer. 
Obviously,  the  two-phase  system  of  a  gas-particle  mixture  reduces  to  a 
single-phase  system  of  a  pure  gas  when  p  =  0.  The  results  for  the  case  of 
P  =  0  should  be  identical  with  the  similarity  solution  for  the  pure-gas 
boundary-layer  equations.  With  p  =  0,  the  numerical  results  for  the 
positions  x  =  1.05,  2.05,  5.05  and  8.05  are  respectively  shown  in  Figs.  18 
to  21.  They  are  compared  with  the  similarity  solution  of  a  pure-gas 
boundary-layer  flow  under  similar  conditions.  Excellent  agreement  between 
these  two  solutions  is  achieved  along  the  whole  plate.  This  comparison 
provides  further  strong  evidence  that  the  present  finite-difference  scheme 
is  correct. 

Once  the  flow  profiles  across  the  boundary  layer  are  known,  the 
boundary-layer  characteri sties  can  be  calculated.  The  three  nondimensional 
quantities  xw,  qw  and  6  as  functions  of  the  distance  x  from  the  leading  edge 
are  shown  in  Fig.  22  to  24.  It  is  found  that  the  curves  for  the  shear 

stress  tw  and  the  heat-transfer  rate  qw  are  nearly  identical  except  in  the 
nonequilibrium  region.  This  is  attributed  to  the  Reynolds  analogy  between 
the  heat-transfer  and  the  shear  stress  [34],  It  is  well  known  that  for  the 
boundary-1  ayer  flow  of  a  pure  gas  on  a  flat  plate,  the  profiles  for  the 
velocity  and  the  temperature  become  completely  analogous  if  the  Prandtl 

number  has  the  value  of  unity.  For  the  dusty-gas  boundary-layer  flows  over 
a  flat  plate,  there  is  no  similar  relation  available.  As  pointed  out 
before,  however,  the  gas-particle  system  behaves  like  a  gaseous  mixture  in 
the  two  limiting  regions:  (1)  In  the  large-slip  limit  near  the  leading 
edge,  the  two-phase  system  responds  as  if  there  were  no  particles.  For  this 
analysis,  the  Prandtl  number  of  the  gas  is  just  assumed  to  be  equal  to  unity 

so  that  the  Reynolds  analogy  holds  in  the  large-slip  region,  or  more 

exactly,  with  zeroth-order  accuracy.  (2)  In  the  small-slip  limit  far 
downstream,  the  two-phase  system  acts  like  a  pure  gas  with  modified 

thermodynamic  properties.  The  modified  Prandtl  number  Pr  is  given  by  [26] 

Pr  =  Pr  i-+  P/g 
1  +  p 

In  the  case  under  consideration  in  this  paper,  a  =  1.0  and  p  =  1.0.  This 

yields  Pr  =  Pr  =  1.0.  Consequently,  the  Reynolds  analogy  also  holds  in  the 

small-slip  region  with  zeroth-order  accuracy.  In  Figs.  22  to  24,  it  is 
interesting  to  note  that  along  every  curve  for  the  characteristic 

quantities,  there  is  an  inflection  point  which  corresponds  to  the  critical 
point.  This  fact  means  that  for  the  boundary-layer  flows  of  a  gas-particle 
mixture,  some  significant  changes  in  the  flow  properties  occur  at  the 
critical  point.  In  fact,  as  mentioned  before,  at  the  critical  point,  the 
two-phase  system  accomplishes  essentially  the  transition  from  the 


nonequi 1 ibrium  flow  to  equilibrium  flow.  In  Figs.  22  to  24,  the 

corresponding  results  for  the  pure  gas  without  particles  are  shown  for 
comparison.  It  is  seen  that  the  shear  stress  and  heat-transfer  rate  at  the 
wall  increase  and  the  displacement  thickness  decreases  owing  to  the  presence 
of  particles.  Owing  to  the  interaction  between  the  particles  and  the  gas  in 
the  two-phase  boundary  layer,  the  tangential  velocity  and  temperature  of  the 
gas  phase  increase  except  at  the  wall  and  at  the  outer  edge,  where  the  same 
boundary  conditions  hold  as  in  the  pure-gas  case.  Compared  with  the 

pure-gas  boundary  layer,  the  velocity  and  temperature  profiles  for  the 
dusty-gas  boundary  layer  have  a  steeper  gradient  at  the  wall  and  a  more  even 
gradient  near  the  outer  edge.  It  is  clear  that,  as  a  result  of  these 

changes  in  the  flow  profiles,  the  shear  stress  and  heat  transfer  increase 
while  the  displacement  thickness  decreases.  The  same  conclusion  was 
obtained  from  the  asymptotic  analysis  [27]. 

Finite-difference  calculations  were  also  done  with  other  types  of 

initial  profiles:  the  zeroth-order  asymptotic  profiles  and  the  extended 
Wu-type  profiles.  In  the  previous  computations  of  this  analysis,  only  the 
first-order  asymptotic  profiles  were  applied.  The  numerical  results  for 
these  three  different  types  of  initial  profiles  are  compared  in  Figs.  25  to 
30.  The  difference  solution  for  the  extended  Wu-type  initial  profiles 
started  at  the  leading  edge  (xQ  =  0.0)  and  the  solutions  for  the  two 
asymptotic  profiles  at  x0  =  0.05.  Figures  25  to  27  give  the  flow  profiles 
of  the  two  phases  at  x  =  0.15,  1.05  and  5.05,  and  Figs.  28  to  30  give  the 
nondimensicnal  boundary-layer  characteristics.  The  results  for  the  extended 
Wu-type  profiles  indicate  that  employing  them  as  the  initial  conditions  can 
result  in  physically  reasonable  solutions  which  agree  very  well  with  the 
solution  for  the  first-order  asymptotic  profiles.  The  results  for  the 
zeroth-order  asymptotic  profiles  show  that  using  the  zeroth-order  asymptotic 
profiles  as  the  initial  conditions  may  cause  some  deviations  in  the  flow 
properties  in  the  quasi -frozen  and  nonequilibrium  regions  but  not  in  the 
quasi-equilibrium  region,  especially  for  the  normal  velocity  of  the 
particles.  The  reason  is  that  the  zeroth-order  asymptotic  profiles  assume  a 
zero  velocity  for  the  normal  particle  velocity  at  some  distance  from  the 
leading  edge  (say,  x0  =  0.05)  and  it  leads  to  some  errors.  With 
increasing  x  along  the  flat-plate,  these  deviations  are  damped  out  and  all 
the  three  initial  profiles  yield  identical  results.  It  should  be  pointed 
out  that  the  possibility  of  using  the  extended  Wu-type  initial  profiles 
leads  to  a  significant  simplification  in  the  computation  procedure,  since  it 
is  not  necessary  to  solve  the  asymptotic  solution  in  the  large-slip  limit. 


8.  CONCLUDING  REMARKS 

The  complete  set  of  nonlinear  partial-differential  equations  for 
compressible,  laminar,  boundary-1 ayer  flows  of  a  dusty  gas  over  a 
semi-infinite  flat  plate  was  solved  using  implicit  finite-difference 
schemes.  The  numerical  solutions  for  the  three  distinct  flow  regimes,  the 
quasi -frozen,  nonequilibrium  and  quasi -equi 1 i bri urn  regimes,  were  obtained 
for  the  case  of  the  Stokes  relation.  The  finite-difference  results  for  the 


two  limiting  cases  of  large  slip  and  small  slip  are  in  good  agreement  with 
the  corresponding  asymptotic  solutions.  The  numerical  examples  indicate 
that  the  present  finite-difference  method  provides  a  useful  technique  for 
studying  two-phase  boundary-layer  flows. 

From  this  analysis,  it  is  shown  that  in  order  to  get  a  finite-difference 
solution  along  the  entire  flat-plate  length  with  the  present  basic 
equations,  it  is  important  to  deal  realistically  with  the  particle  density 
after  the  critical  point.  The  assumption  that  the  gas-particle  system  is 
treated  as  a  binary  gas  with  the  given  mass  ratio  of  the  components  ( p) 
after  the  critical  point  represents  a  practical  approach  for  the  case  where 
diffusion  is  the  main  controlling  process  in  the  region  x  >  xcr^ .  It  yields 
physically  reasonable  results. 

For  the  x-momentum  and  energy  equations  of  the  particles,  both 
four-point  and  six-point  schemes  can  be  used.  The  numerical  computations 
indicate  that  the  six-point  scheme  leads  to  better  results  especially  in 
the  nonequilibrium-flow  region,  since  it  has  an  accuracy  of  second  order. 
When  using  the  six-point  scheme  it  is  necessary  to  obtain  the  compatibility 
conditions  as  additional  boundary  conditions.  Fortunately,  the 
compatibility  conditions  for  the  tangential  velocity  and  temperature  of  the 
particles  are  very  simple  after  the  critical  point. 

A  comparative  study  of  the  three  different  initial  profiles  (the 
first-order  asymptotic,  zeroth-order  asymptotic  and  extended  Wu-type) 
indicates  that  all  three  types  of  initial  profiles  can  be  used  for  the 
finite-difference  solution  but  the  zeroth-order  asymptotic  profiles  might 
cause  some  errors  in  the  near  leading-edge  and  transition  regions.  It  is 
suggested  that  the  extended  Wu-type  initial  profiles  can  be  used,  since  they 
lead  to  a  significant  simplification  in  the  numerical  procedure. 

The  numerical  results  presented  in  this  report  indicate  that  the  effects 
of  the  particles  on  the  boundary-layer  flows  are  considerable  and  that  the 
modification  in  the  flow  properties  owing  to  the  presence  of  particles 
includes  an  alteration  of  the  flow  profiles,  an  increase  in  the  shear  stress 
and  heat  transfer  at  the  wall  and  a  decrease  in  the  displacement  thickness. 
These  are  the  major  features  found  in  compressible,  laminar  boundary-layer 
flows  of  a  gas-particle  mixture. 
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FIG.  4  FLOW  PROFILES  IN  LARGE-SLIP  REGION  (x  =  0.055). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

O  DIFFERENCE  SOLUTION  FOR  GAS;  O  DIFFERENCE  SOLUTION  FOR  PARTICLES; 
---  ASYMPTOTIC  SOLUTION  FOR  GAS;  ASYMPTOTIC  SOLUTION  FOR  PARTICLES. 
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FIG.  b  FLOW  PROFILES  IN  LARGE-SLIP  REGION  (x  =  0.105). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

O  DIFFERENCE  SOLUTION  FOR  GAS;  O  DIFFERENCE  SOLUTION  FOR  PARTICLES; 
---  ASYMPTOTIC  SOLUTION  FDR  GAS;  ASYMPTOTIC  SOLUTION  FDR  PARTICLES. 
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FIG.  6  FLOW  PROFILES  IN  LARGE-SLIP  REGION  (x  =  0.305). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

°  SOLUTION  F0R  GAS;  O  DIFFERENCE  SOLUTION  FOR  PARTICLES; 

---  ASYMPTOTIC  SOLUTION  FOR  GAS;  ASYMPTOTIC  SOLUTION  FOR  PARTICLES. 
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FIG.  8  FLOW  PROFILES  IN  MODERATE-SLIP  REGION  (x  =  1.05). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE 
O  GAS;  O  PARTICLES. 
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FIG.  9  FLOW  PROFILES  IN  MODERATE-SLIP  REGION  (x  =  2.05). 

(a)  TANGENTIAL  VELOCITY;  { b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 
O  GAS;  O  PARTICLES. 


/ivjV’.V;  - 


FIG.  1 


7  COMPARISON  BETWEEN  FLOW  PROFILES  RESULTING  FROM  FOUR-POINT  AND  SIX-POINT 
SCHEMES  IN  SMALL-SLIP  REGION  (x  =  10.05). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 


FOUR-POINT  SCHEME:  O  GAS;  O  PARTICLE; 
SIX-POINT  SCHEME:  ---  GAS; - PARTICLE. 
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FIG.  18  COMPARISON  BETWEEN  FLOW  PROFILES  RESULTING  FROM  FINITE -DIFFERENCE 
SOLUTION  WITH  0*0  AND  SIMILARITY  SOLUTION  AT  x  =  1.05. 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

O  FINITE-DIFFERENCE  SOLUTION  WITHOUT  PARTICLES;  —  PURE-GAS  SIMILARITY 
SOLUTION. 
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FIG.  20  COMPARISON  BETWEEN  FLOW  PROFILES  RESULTING  FROM  FINITE-DIFFERENCE 
SOLUTION  WITH  0  =  0  AND  SIMILARITY  SOLUTION  AT  x  =  5.05. 


(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 


O  FINITE-DIFFERENCE  SOLUTION  WITHOUT  PARTICLES;  -  PURE-GAS  SIMILARITY 

SOLUTION. 
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FIG.  21  COMPARISON  BETWEEN  FLOW  PROFILES  RESULTING  FROM  FINITE-DIFFERENCE 
SOLUTION  WITH  0  =  0  AND  SIMILARITY  SOLUTION  AT  x  =  8.05. 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

O  FINITE-DIFFERENCE  SOLUTION  WITHOUT  PARTICLES;  -  PURE-GAS  SIMILARITY 

SOLUTION, 
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FIG.  26  FLOW  PROFILES  WITH  THREE.  DIFFERENT  TYPES  OF  INITIAL  PROFILES  IN 
NONEQUILIBRIUM  FLOW  REGION  (x  =  1.05). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

ZERQTH-ORDER  ASYMPTOTIC  PROFILES:  O  GAS;  O  PARTICLES; 

EXTENDED  WU-TYPE  PROFILES:  A  GAS;  A  PARTICLES; 

FIRST-OROER  ASYMPTOTIC  PROFILES:  -  —  GAS;  —  PARTICLES. 


FIG.  27  FLOW  PROFILES  WITH  THREE  DIFFERENT  TYPES  OF  INITIAL  PROFILES  IN 
NEAR-EQUILIBRIUM  FLOW  REGION  (x  =  5.05). 

(a)  TANGENTIAL  VELOCITY;  (b)  NORMAL  VELOCITY;  (c)  TEMPERATURE. 

ZEROTH-ORDER  ASYMPTOTIC  PROFILES:  O  GAS;  O  PARTICLES; 

EXTENDED  WU-TYPE  PROFILES:  A  GAS;  £  PARTICLES; 

FIRST-ORDER  ASYMPTOTIC  PROFILES:  —  GAS; —PARTICLES. 
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APPENDIX  A 


DERIVATION  OF  THE  FINITE-DIFFERENCE  EQUATIONS  WITH  A  SIX-POINT  SCHEME 


The  momentum  equation  of  the  gas  phase  is 


pu  &u  +  py  =  n  +  4#  ili  +  pn(un  -  u)pD 


dx 


*  '  ^  dT  ay  1*  KPVUP 


or 


pu  8y.  +  (pv  -  H*T y)  -  PpUpPD  -  PpUpD 

dx  °y  ay  2  y  y 


( A- 1 ) 

(A-2) 


where 


I**  =  IP, 
dT 


T  aT  =  Tn+1  +  (K2-1)Tn  ~  k2tp-1 

'  i  •  ~~  ~ 


y  ay  (1+K)  Ayn 


Substituting  the  expressions  for  difference  quotients  (3.1)-(3.4),  Eq. 
(A-2)  becomes 

(P^m+e.n  ^  (um+l,n  "  um,n)  +  < *  “  ^'Vm+e.n  {(1+^  Ay”  Kn+l.n+1 


+  (K2_l)um+l,n  '  KVl,n.ll  +  Vn  > 


Pm+e,n  ^ — - -  Kn+l.n+l  “  (K+1)um+l,n  +  Kum+l,n-J 

(l+K)Ayn2 


+  2(1-8)K  A2g  l 

—  _  u  um,n 1 

(1+K)  Ayn  2 

=  (PpupP°)m+9,n  "  [0(  PpP°)m+l,n  unH-l,n  +  (1"e)(  PpPD)m,n  um,n^ 


(A-3) 


Multiplying  by  Ax  and  rearranging  Eq.  (A-3),  the  finite-difference  form  of 
the  momentum  equation  can  be  obtained: 


An  um+l ,n+l  +  Bn  um+l,n  +  Cn  Jm+l,n-l  =  Dn 


(A-4) 


where  the  coefficients  have  the  following  expressions: 


An  =  an^pv  “  p  ^y  Ws,n  "  cn  Pm+6,n 


Bn  =  (P^m+G.n  +  an(K2-l)(pv  -  ^Ty)m+0fn  +  cn(K+l)  i^e.r,  +  9Ax(  f^ii))m+1  >n 


^n  =  “anK2(pv  "  P  ^yWs.n  "  cnK|Jm+0,n 


(A-5) 


Bn  =  [(Pu)m+9,n  “  (I- ®)  Ax(  P^)m,n  -lum,n  “  ^(pv  "  p  ^y^m+O.n^m.n 


+  ^nPm+G.n  A2um,n  +  Ax(  ppupPD)m+0,n 


The  energy  equation  of  the  gas  phase  is 


+  pv  —  =  -H-  +  J_  lii  (i*L}2  +  Ec  p(—  )2  +  Ec  pu  [(uD  -  u) 

9x  9y  Pr  Qy  2  Pr  dT  'ay  'ay  ^  p 


+  (VP  '  V)2]^  +  3?F  PP(TP  "  ™U 


( A— 6 ) 


PU  +  (pv  -  T y)  -  £-  —  =  Ec  v&L)2  +  Ec  Pp[(up  -  u)  2 

ax  Pr  y  ay  Pr  av2  ay  v  v 


+  1  (vp  -  v)2]mD  +  PpTppNu  -  A-  PpTpNu 

00 


(A-7) 


Using  a  similar  procedure,  Eq.  (A-7)  is  discretized  as 


( pu  We,n  ^  ( rm+l  ,n 


Vn)  +  (pv  '  ff  Jy  Ue,n^(l+K)Ayn  ^m+1’n+1 


Conti nued 


ii»««  wi  mranmnrwnvRvwv  w 


TW 


+  (K2-l)Tm+l,n  ‘  K2Tm+l,r>-lJ  +  ATm,n^ 


-  (#-).„{- 


2  6K 


Pr  m+0,n  (i+K)Ayn2 


t^m+l,n+l  "  ^+^^m+l,n  +  ^m+l,n-l  ] 


+  ^Lb9)K  ^  } 

,  ,  „  m,n  > 

(1+K)  Ayn2 

{Ec  p  uy2  +  Ec  Pp [ ( u p  “  u)  2  +  -i-  (vp  -  v)  2]mD  +  -1-  f^TpMNu  }m+0>n 


Wl,n^im-l,n  +  ^-0)  Pp^  )m,nTm,n  ]  (A“8) 


3Pr 


The  difference  equation  for  the  energy  conservation  of  the  gas  is 


2  2  2  2 
An  Tm+l,n+l  +  Bn  ^m+l,n  +  ^n  ^m+l,n-l  ~  Dn 


(A-9) 


where 


An  =  an^pv  “  p~  ^yJm+O.n  "  cn(p“l™+0,n 


Pr 


Bn  =  (P^m+O.n  +  an(K2-l)(pv  -  Ty)^^  +  cn(W)^Ue,n 


Pr 


+  9AX^3?F 


(A- 10) 


^n  =  -an^2(pv  “  p“  ^y)m+0,n  “  cn^prli+0,n 


Pr 


Dn  t(  P'j)m+e>n  ■  (1-0)  Ax  (— L-  Pp  pNu  )m>n  ]Tmj  n  -  bn(pv  -  j~-  Ty  ATm^n 

+  dn(#rUo,n  A2Tm,n  +  **{Ec  puy2  +  EcPp[up  -  u)2  +  -i_  (vp  -  v)2]pD 


Continued 


+  o5T  pp  p^u  ta+e.n 


The  x -momentum  equation  of  the  particle  phase  is 


UP  ^  +  VP  ^  =  "(UP  "  U)^  =  UM°  ‘ 


( A- 11) 


Discretizing  Eq.  ( A- 11)  yields 


UPm+0,n  Ax  ^  Pm+l,n  Pm,n  ^  Pm+  9,n  ^(l+K)  Ayn  ^  Pm+l,n+l 


Tisfc- 


(uP°)ni+9,n  '  [^PDJm+l.nX+i.n  +  (1-9)(  tjD)m,nUpm#n  ]  (A'12) 


Therefore 


where 


A"  UPm+l  ,n+l  +  Bf1  UPm+l  ,n  +  C"  UPm+l,n-l  =  °n 


An  =  a"  VPm+0,n 


*"  =  V.,»  +  a"(KM)Ve.r,  +  ^(MD)-+l.n 


( A- 13) 


( A-14) 


Cn  ::  'anK2vpm+e>n 

°"  =  [uPm+9,n  "  (1-e)Ax(M0)-.n]X,n  '  b"VPm+9,n  Xn  + 


The  energy  equation  of  the  particle  phase  is 


WyCNC 


*  s. 


( A- 15 ) 


ST  5T 

un  — 1  +  vn  —E  =  -  .JL  (T  -  T)pNu  =  -SL  T uNu  -  _JL  Tn(iNu 

P  9x  P  3y  3Pr  P  3Pr  3Pr  P 


Then 


Pm+0,n  Ax  ^  Pm+l,n  Pm,n^  VPm+6,n  ^(1+K)  Ayn  ^Pm+l,n+l  ^  ^Pm+l,n 


-  K2Tn  ]  +  -All _  AT-.  } 

Pm+l,n-l  (l+K)Ayn  Pm,n 1 


’  (WF  T“NU)™^."  '  [6(W  *  TPmtl,n  +  (1-9)(3?F  **  W»  TP„,J 

(A-16) 

Consequently, 


B  8  8 

A  T  +  B_  T  +•  C_  T_ 

n  Pm+l,n+l  n  Pm+l,n  n  Pm+l,n-l 


=  0. 


(A-17) 


where 


An  =  anVPm+0,n 


Bn  =  u0  +  an(K?  1 )  vn  +  0Ax  (— —  pNu  L.  i  n 

n  Pm+0,n  n  Pm+0,n  3Pr  ^+1‘n 


C  =  -anK2  v 

n  n  Pm+0,n 


( A— 18) 


°n  =  tuPm+0,n  "  (1'9)Ax(3FF  ^‘"^Pm.n  "  b"VPm+0,n  ATPm,n 


+  T^U0,n 


The  continuity  equation  of  the  particle  phase  can  be  vritten  in  the 
different  forms  as  follows: 


Un  -ZR  +  +  pn(Jt_l)  =  o 

P  dx  P  b/  P  5x  dy 


( A— 19) 


Up  f  (— )  *  vp  »  (JL)  -  (L)A  *  £i>)  =  0 
P  dx  P  dy  Pp  Pp  dy 


The  third  form  can  be  discretized  as 


,  1_  rrL.)  .  (!_]  1  +  v  {  6 _  [(!_) 

P(TH-0,n  Ax  pp  m+l,n  Pp  m,n  Pm+0,n  (l+K)Ayn  Pp  m+l,n+l 


+  <K2'1)(i;),»*l,n  •  +  (T^ 


-  [9(-i-)m*l,n  +  (uPm+1>n  '  uPm>„'  +  (l*K)4yn"  4vPnH-l,n 


+  1- :J. _ Avn  }  =  0 

(1+K)  Ayn  pm,n 


( A-20) 


Then  the  difference  equation  becomes 


M-^Wl.n+i  +  M*  L1>n  +  Cn^p  Ul,n-1 


(A-21] 


The  coefficients  in  Eq.  (A-21)  are: 


Af1  =  anVPm+9,n 


B  =  un  +  an (K 2-l ) vn  -  9(un  -  uD  )  -  a  Av„ 

n  pm+9,n  n  Pm+0,n  Pm+l,n  Pm,n  n  Pm+0,n 


Cn  =  -anK  2v 
n  n  Pm+0,n 


(A-22; 


*  w"*  »*■  >.*»■ ,  ’» k.’ 

“fk/h v"V ■  ..•r . if -1 « r .  1  .■fV’V' , 


APPENDIX  B 


DERIVATION  OF  THE  FINITE-DIFFERENCE  EQUATIONS 


'I 

4 


WITH  A  FOUR-POINT  SCHEME 


The  continuity  equation  of  the  gas  phase  is 


* 

$ 

‘ft 


•if1 


iL  pu  +  JL  pv  =  o 
ax  ay 


(B-l) 


Using  the  expressions  for  the  finite-difference  quotients  (3.5)-(3.7),  Eq, 
(B-l)  becomes 


[(Pu)m+l,n  "  (P^m.n  +  i^W^n-l  "  (P^m.n-lJ 


a 


f 


1 


Vi 

I 


A*n-1 


[(Pv)m+i  n  “  ( pv)nH-l,n-l  1  +  -  ^.(pv)m  n  -  (pv)m,n-l]  =  ® 

A*n-1  (B-2) 


Therefore 


(pv)m+l,n  =  (Pv)m+l,n-l  "  ^P^m.n  “  (pv)m,n-l3 
”  P0a7~  "  ^pu^m,n  +  ^Wl.n-l  “  (P^m.n-ll  (B-3) 


The  x-momentum  equation  of  the  particle  phase  is 


t* 

V 

$ 

,y 


I 

A 

-«* 

A 


aun  auD  .  .  _ 

UD  +  Vn  ~T^  =  ” ^ UD  “  U)P° 

p  ax  p  ay  p 


With  the  quotient  expressions  (3.8)-(3.10),  Eq.  (B-4)  becomes 


i~  (u 


-  un  )  +  v 


14-  <», 


(B-4) 


Pm+9,n  Ax  Pm+l,n  Pm,n  Pm+0,n  Ayn  Pm+l,n+l  pm+l,n 


Continued 


‘v.v 

ViY 

•1$ 


I  V, 

d 

.'fV 

Ati*i 


+  ( 1  -.0)  (up  -  Up  )} 

Ayn  pm,n+l  pm,n 

=  -[<XM»Wl.n  Vl,»  *  (1-e)(|i,>™,n  V"1  * 


Therefore,  the  difference  equation  can  be  written  as 


(B-5) 


K\V 

‘vHS 


■*pVi 

I 


where 


A"  Upm+1 ,n+l  +  Bn  Upm+1 ,n  =  C" 


A3  =  V 

n  Ayn  pm+  9,n 


p  =  u  -  -  -  v_  +  9Ax(pD)m+i  n 

"n  pm+9,n  Ayn  pm+8,n  * 


(B-6) 


(B-7) 


i 

fs 


v.tb 


i 

is; 


7M 

1 


I 


c3  .  -[(1-9)  Ax  v  ]u  +  [u_  +  ill V  „ 

n  L  Ayn  pm+9,nJ  pm,n+l  pm+9,n  Ayn  pm+0,n 

-  (l-0)Ax(MO)mtr,)jpm  +  Ax(uMO)m+0,n 


The  y-mocnentum  equation  of  the  particle  phase  is 


“p  ^ *  vp  I2  =  -(vp  ' v)'10 


Substituting  expressions  (3.8),  (3.9)  and  (3.11)  into  Eq.  (B-8), 


^  r  0  /  ' 

Upm+9,n  A^  (Vpm+l ,n  '  Vpm,n)  +  Vpm+9,n  ^y^'  (Vpm+l,n  Vpm+I,n-1‘ 


(B-8) 


+  ..lli-  (vD  -  vp  )} 
Ayn_1  pm,n  pm,n-l 


m 


•VS 


-  V1>n  *  vpm<n i  * 


(B-9) 


Then,  the  finite-difference  form  for  the  y-momentum  equation  becomes 


A"  VPm+l,n  +  Bn  VPm+l,n-l  =  °n 


"  =UPm+e,n+^TVPm+e,n  +  9Ax(^Wl,n 


( B— 10) 


B1*  =  -  v 

n  Ayn-1  pnH-8,n 


(B-ll) 


c"  '  ^UPtn+0,n  '  VPm+9,n  '  (1'6)  MD)m.n 

*  [i^T  Ve.JVn-l  *  4XlV ^ "+e>" 


The  energy  equation  of  the  particle  phase  is 


uD  HE  +  v-  5e  =  -  _5L  (T  -  T)pNu 
P  9x  P  dy  3Pr  P 


(B— 12) 


Using  a  procedure  similar  to  the  x-momentum  equation,  the  finite-difference 
equation  for  particle  energy  conservation  can  be  derived  as  follows: 


i-(T 


-  V  >  ♦*«.. 


-  Tr 


Pm+0,n  Ax  Pm+l,n  Pm,n  Pm+9,n  Ayn  Pm+l,n+l  Pm+l,n 


) 


+  ^(TPm,n+l-TPm,n)} 

**W.n  Vl.n  +  (1-9)(3F  TPm,n]  +  <*T*“U.n 

(B— 13) 


Af1  TPm+l ,n+l  +  Bn  TPm+l,n  =  Cn 


( B— 14) 


8,3 


A 5  =  v 
n  Ayn  pm+9,n 


b!  =  u 


n  Pm+9,n  Ay 


"  VPm+0,„  *  eSx<W’J(u)”+‘." 


c5  =  v  It  +  ru  +  LkiL^  v 

n  Ayn  Pm+9,n  Pm,n+1  '  Pm+9,n  Ayn  Pm+9,n 


The  continuity  equation  of  the  particle  phase  is 


—  a,u_  +  —  p_  v,,  =  0 
ft*  HP  P  ftw  HPVP 


(B-16 


aPn  aPn  .  9Un  3vn  . 

u  — P  +  v_  — E-  +  p_( — P-  +  — P-)  =  0 
p  ax  P  ay  Pvax  ay  ' 


Usi^g  the  quotient  expressions  (3.8)-(3.10) ,  Eq.  ( B- 17)  becomes 


4:  Cp, 


Pn_  J  +  VD_.  A  (pD 


Pm+9,n  Ax  Pm+l,n  Pm,n  J  Pm+9,nLAyn  Pm+l,n+l  Pm+l,n 


+  (pPm,n+l  "  PPm,n)]  +  [9pPm+l,n  *  (l‘9)Vn^  (uPm+l,n  "  X.J 


+  —  (v  -  V  )  +  ^  (v_  -  V  11=0 

Ayn  Pm+l,n+l  Pm+l.n  Ayn  Pm,n+1  pm,n 


Then,  the  finite-difference  equation  can  be  written  in  the  following  form 


A°  PPm+l  ,n+l  +  Bn  PPm+l  ,n  =  C" 


APPENDIX  C 


DERIVATION  OF  THE  RELATIONS  FOR  SHEAR  STRESS,  HEAT  TRANSFER 
AND  DISPLACEMENT  THICKNESS 


Characteristic  quantities  of  boundary-layer  flows  are  defined  as 


is  =  k:  (£J) 

w  Vt  Qy*  w 

(c-1) 

4v5  = 

^  w  by*  w 

(C-2) 

6*  -  /*  (l  -  J^)dy* 

Q  U 

Q  r  CD  OO 

(C-3) 

Correspondingly,  the  nondimensional  parameters  take  the  following  form: 


*w 


— vUe* 
Pi  u*2 


where 


6  = 


vUe 

K 


oo 


(C-4) 

(C-5) 

(C-6) 


From  the  definition  of  the  nondimensional  parameters,  Eq.  (2.11),  the 
derivatives  of  the  gas  velocity  u*  and  temperature  T*  can  be  given  as 


and  the  integration  (C-3)  becomes 


/  (1  -  )dy*  =  !  (1  -  pu)dy  (C-9) 

o  P®  u»  yRe„,  o 


Therefore,  the  nondimensi onal  characteristics  of  the  boundary-layer  flows 
can  be  expressed  as. 


\  =  Pw  ) 
w  w  5y  w 


(C-10) 


Pw  ( 3T  j 
Pr  Ec  ^by  w 


(C-ll) 


where 


6  =  /  (1  -  pu)dy 
o 


i*  2 


and 


Ec  = 


c*  T* 

1  ot> 


(C-12) 


In  order  to  calculate  the  derivatives  at  the  wall  with  the 
finite-difference  solutions,  the  gas  velocity  u  near  the  wall  can  be 
approximated  by  a  cubic  polynomial: 


u  =  au  +  buy  +  cuy2  +  duy3 


( C-13) 


At  the  four  grid  points  nearest  the  wall,  the  values  of  gas  velocity  are 
known  and  are  equal  to  u^ ,  U2,  U3  and  u^,  respecti vely.  Then,  there  are 
four  equations  which  can  be  used  to  obtain  the  coefficients  au,  bu,  cu  and 


y  1  =  0: 


U1  =  au 


( C- 14) 


y  =  Ay  : 

2  1 


U2  =  au  +  buAyi  +  cu  Ay ! 2  +  d^3 


(C-15) 


y3  =  (K+l)Ayi:  u3  =  au  +  b^K+l^  +  cu(K+l)2Ay12 


+  du(K+l)3Ayi3 


(C-16) 


y4  =  (K  2+K+l )  Ay^ :  U4  =  au  +  b^K^K+l)  Ayj  +  cu(K  ^K+l)  2Ayj  2 


+  du(K2+K+l)  3Ayt  3 


Solving  the  system  of  simultaneous  equations  (C-14) 
elimination  method,  the  value  of  coefficient  bu  is  obtained: 


(C-17) 


(C-17)  by  an 


bu  -  [(u2  -  uo  -  +  Alui..] 

K2Ayi  *  K(K+1)  K (K  2+K+i ) 2 


(C-18) 


For  the  gas  temperature  T,  there  is  a  similar  relation  for  coefficient 


bT  ■  [(T,  -  T.)  -  I|_:  .ll  +  .Tj,  ;_T.l  ]  (C-19) 

K24yj  K(K+1)  K(K  2+K+l) 2 


Since  the  gas  velocity  at  the  wall  vanishes,  uj  =  0  in  the  expression 
(C-18),  therefore,  the  nondimensional  shear  stress  at  the  wall  can  be 
determined: 


V  -  [u2  -  _ ] 

K2Ay1  K(K+1)  K(K2+K+1)2 


(C-20) 


and  the  corresponding  heat  transfer  is 


•  __  _  Jw_  (K^KlU  [(T  .  T  )  .  1.3  ;  T.l  +  __4_2._1 —  j  {c_21) 

Pr  Ec  K2Ay1  K(K+1)  K(K  2+K+l) 2 


As  for  the  three-point  difference  formula,  the  integration  I  is 
calculated  by 


I  =  /  F(y)dy 
o 


(C-22) 


Within  every  small  integration  region,  the  integrated  function  F(y)  can  be 
approximated  by  a  quadratic  parabola  p(y)  through  three  consecutive  grid 
points  (m+1,  n-1),  (m+1,  n)  and  (m+1,  n+1): 


p  (y )  =  An(Llln)2  +  Bn(LlZn)  +  cn  (c-23) 

n  Ayn  n  Ay^  n 


The  quadratic  polynomial  p(y)  satisfies  the  conditions: 


P(yn_l)  1 , n - 1 »  "  Fnt+l,n»  P^n+l^  Fm+l,n+l  ( C-24) 

Then  a  system  of  simultaneous  algebraic  equations  for  the  coefficients  An, 
Bn  and  Cn  can  be  constructed: 

An  "  KBn  +  K^n  =  K2Fm+l,n-l»  ^n  =  ^m+l.n’  An  +  Bn  +  Cn  Fm+ijn+i 

(C-25) 


Consequently,  the  solutions  of  the  above  equations  can  be  found  as 


a  _  K  p  L/p  .  K 2  p 

nn  r m+1, n+1  "  Nrm+l,n  jyyy  r rn+ 1 , n - 1 


“f  p m+1, n+1  +  ^“^m+l  ,n  "  ^m+l.n-l 


K+l 


(C-26) 


pm+l,n 
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Substituting  the  coefficient  expression  (C-26)  into  Eq.  (C-23),  replacing 
the  integrated  function  F(y)  by  p(y),  and  integrating  Eq.  (C-22)  in  every 
small  interval  Ayn-1  and  making  the  summation,  the  three-point  difference 
formula  for  nonequal  intervals  is  derived: 


l.’l  t^-An  -4SSM-»n*  ^-iCj 

2  3K2  2K 


.  V  4yn-i  r3K+2  F  ,  3K+1  F  _  1  c  .  .  1  (r.?n 

-  I  t-m+i.n-i  +  nr-  m+1»n  kTk+iT  m+l.n+l J  ; 


k(k+i; 


APPENDIX  D 


COMPUTER  PROGRAM  FDBLEP 


The  program  FOBLEP  for  solving  finite-difference  boundary-1 ayer 
equations  over  a  flat  plate  for  a  dusty  gas  is  written  on  the  Perkin-Elmer 
3250  system  at  UTIAS.  The  main  notations  used  in  the  program  are  listed  and 
explained  below. 


U(I)  - 
V(I)  = 
T  ( I )  = 
R0(I)  = 

UP(I)  = 
VP(I)  = 
TP ( I )  = 

ROP(I)  = 
SHEAR  = 
HEAT  = 


u  , 

at 

grid 

poi  nt 

(m+1. 

n) 

V  . 

at 

grid 

poi  nt 

(m+1. 

n) 

T 

at 

grid 

point 

(m+1. 

n) 

P  ■ 

at 

grid 

point 

(m+1. 

n) 

UP 

at 

grid 

point 

(m+1. 

n) 

VP 

at 

grid 

point 

(m+1, 

n) 

TP 

at 

grid 

point 

(m+1, 

n) 

PP 

at 

grid 

point 

(m+1. 

n) 

at 

grid 

line 

(m+1) 

% 

at 

grid 

line 

(m+1) 

6  at  grid  line  (m+I) 


THICK 


vl'l 

w1 


m 


i 

i?  **> 


maximum  value  of  n 


initial  value  of  x 


maximum  value  of  x 


critical  value  of  x 


non 


- - - - - - 

C  THIS  PROGRAM  SOLVES  THE  GAS  AND  PARTICLE  PARAMETERS  FOR 

C  THE  BOUNDARY  LAYER  EQUATIONS  OF  A  DUSTY  GAS  OVER  A  SEMI- 

C  INFINITE  FLAT  PLATE  BY  MEANS  OF  FINITE  DIFFERENCE  METHOD 

C  AND  GIVES  CHARACTERISTIC  QUANTITIES  OF  BOUNDARY  LAYER  FLOW 

C  (IMPLICIT  SCHEMES  FOR  TWO  PHASES.  NONITERATION  PROCEDURE) 

c  -  -  -  —  —  —  —  -  -  -  -  -  - - —  -  -  -  - - 

c  -  -  -  —  -  -  -  -  —  -  -  -  -  -  -  - - -  -  -  - - 

C  MAIN  PROGRAM 

C  - - - - - - - - - - - - - —  -  -  — 

- - - - 

IMPLICIT  REAL*B(A-H.  O-Z ) 

REAL*8  KE.  KE1 ,  KE2,  KE3,  KE4,  MUI ,  MDR.  MNU,  MUW 
DIMENSION  A(  100).  B  (  100  ) ,  C  (  100) ,  D<  100) 

DIMENSION  U<  100).  V<  100  ) ,  T (  1 00 ) ,  RO <  1 00 ) 

D I  MENS  I  ON  UP(100>,  VP(IOO),  TP(IOO). ROP (  100) 

DIMENSION  RUB < 100)  ,  RVB < 100 ) 

DIMENSION  RU<  100),  CR<  100).  CRB  <  100) 

DIMENSION  DMDf 100)  ,  DMN < 100  > . DMDB ( 100) ,  DMNB < 100 ) 

DIMENSION  UD(  100),  VD<  100),  TD(  100) 

COMMON  /G/  PR.  WN,  EC 
COMMON  /HI/  HE. K£1 . KE2, K£3 
COMMON  /H2/  CA.CB.CC,  CD 

C  -  -  - - - - - - - — ■  — 


1020 

FORMAT 

(3E13  6) 

1030 

FORMAT 

(5E13  6) 

1050 

FORMAT 

(IX, 3E13  6) 

1060 

FORMAT 

(IX, 5E1 3  6) 

1070 

FORMAT 

(21  5,  2E13  6  ) 

1080 

FORMAT 

( IE  16  8, 7E1 5  8) 

1090 

FORMAT 

(IX,  215,  2E13  6) 

1 100 

FORMAT 

(IX, 8E15  8) 

1110 

FORMAT 

( IX,  215,  1E13.  6,  3E15  8) 

SET  THE  BASIC  PARAMETERS 


READ  (7,1020)  UW,  TW,  BETA 

READ  (7,1020)  EC,  PR,  WN 

READ  (7,1030)  DX.  D Y,  EPS,  KE,  CITA 

WRITE  (B,  1050)  UW,  TW,  BETA 

WRITE  (6,1050)  EC,  PR,  WN 

WRITE  (8,1060)  DX,  DY,  EPS,  KE,  CITA 

READ  (7,1070)  MEND,  N,  XSTA,  XMAX 

READ  (7,1080)  ( U<  I  ) ,  V  (  I  ) ,  T  (  I ) ,  R0(  I  ) , 

1  UP(  I  ),  VP(I  ) ,  TP<  I  >.  ROP<  I  >,  1  =  1, N) 

WRITE  (8,1090)  MEND,  N,  XSTA,  XMAX 

WRITE  (8,1100)  <U<  I  ),  V<  I  >,  TCI ) ,  R0<  I  ),  UP(I  ) ,  VP  (  I  ).  TP(I  ) ,  ROP  (  I ) 
1  1=1 , 10) 


0.3 


n  r>  r>  non 


SET  THE  NUMERICAL  PARAMETERS 


KE1=KE  +  1.  D+OO 

KE2=KE#KE 

KE3=KE2-1  D+00 

KE4=KE 1+KE2 

C I TA1  =  1 .  D+OO-CITA 

C I TA2=C ITA 1 /C I TA 

C I TA3=2.  D+00*C ITA 

C1TA4=CITA3-1  D+OO 

C I  TA5=2.  D+00*C ITA1 

C I  1 =C I TA*D X 

C 1 2=C  I TA1 *DX 

Cl  3=2  D+OO *C II 

Cl 4=2.  D+OO *C I 2 

CK1=KE4/KE2 

CK2=KE*KE1 

CK3=KE*KE4*KE4 

CK4=( 2  D+00+3  D+00*KE)/KE1 

CK5=< 1 .  D+00+3  D+00*KE)/KE 

CK6=1 .  D+00/CK2 


CALCULATE  THE  FLOW  PROFILES  AT  THE  NEXT  GRID  LINE  (M+l ) 


MUU=TW**WN 
XCRI=1 .  D+OO/MUW 
X=XSTA 

10  CONTINUE 

DO  950  L=1 ,  10 

DO  900  M= 1 ,  MEND 

N1=N-1 

N2=N-2 

X=X+DX 

IF  (X  GT.  XMAX)  GO  TO  999 


C  STORE  THE  VALUES  OF  GAS  VELOCITY  AT  THE  PREVIOUS  LINE  (W 


DO  20  1=1, N 
20  RU ( I )=U< I > 

C  -  -  - - - - - - - - - - - - - - 

C  SOLVE  THE  NEW  GAS  TANGENTIAL  VELOCITY  U  AT  THE  NEXT 

C  GRID  LINE  < M+ 1 )  USING  SIX-POINT  IMPLICIT  SCHEME 

C  -  -  - - - - - - - - - - - - - - 

C  SET  THE  BOUNDARY  CONDITIONS  AT  THE  LINE  <M+1) 

C  -  -  - - - - - - - - - - - - - - 

U<  1  >=UW 
U(N)  =  1 .  OD+OO 


n  n  o  *-  non  non  oooo 


SET  THE  FLOW  PARAMETERS  AT  THE  GRID  POINT  <N*1  )  TO  CALCULATE 
THE  FIRST  AND  SECOND  ORDER  DIFFERENCES  AT  Tt£  GRID  LINE  (M) 


U(N+1)-U(N> 

T(N+1)-T<N) 


CALCULATE  THE  INITIAL  VALUE  OF  SOME  COEFFICIENTS 


CALL  VALU1 (CIl, CI2,  CI3,  CI4,  DY,  DY1 ) 


CALCULATE  THE  COEFFICIENT  MATRIX  ELEMENTS 


CALL  PARA1  (  1,  U,  V,  RO,  RUI,  RVI  ) 

RUB ( 1 ) =RUI 
RVB ( 1 >=RVI 

CALL  PARA2<  1,  T,  MUI ,  DMUI,  MDR,  MNU) 

CALL  PARA 5  ( 1,  DX,  MDR,  MNU,  DMDB,  DMNB,  U,  V,  T,  UD,  VD,  TD) 

DO  110  1=2, N 

CALL  PARA1  <  I,  U,  V,  RO,  RUI,  RVI  > 

RUB  < I ) =RUI 
RVB  < I ) =RVI 

CALL  PARA3  <  I,  KE2,  KE3,  U,  DEUI ) 

CALL  PARA3  ( I,  KE2,  KE3,  T ,  DET I  ) 

DY1=DY1*KE 
DT I=DETI/DY1 

CALL  PARA4  (  I,  KE,  KE1,  U,  DDEUI  ) 

CALL  PARA2<  I,  T,  MUI,  DMUI.  MDR,  MNU) 

CALL  P  ARA5  <  I,  DX,  MDR,  MNU,  DMDB-  DMNB,  U,  V,  T,  UD,  VD,  TD) 

ROMD=ROP( r  >*DMDB< I ) 

UDI=ROMD*UP(I) 

CR  <  I  )=C ITA*RQMD 
CRB  < I ) =CITA1*R0MD 
DMUI=DMUI*DTI 
CALL  VALU2  <KE,  KE2 ) 

10  CALL  C0EF1  (  I,  U.  RUI,  RVI,  MUI,  DMUI,  DEUI,  DDEUI,  UDI,  CR,.CRB,  A,  B,  C,  D> 
CALL  C0EF2  (U,  A,  B,  C,  D) 


SOLVE  ThC  GAS  TANGENTIAL  VELOCITY  PROFILE  BY  THOMAS  ALGORI Trtt 


DO  120  1=2.  N1 

120  CALL  TH0M1  (  I,  A,  B,  C,  D) 

DO  130  1=1 , N2 
J=N— I 

130  CALL  TH0M2  (  J,  U,  A,  B,  D) 

C  -  -  - - - - ______ - - - - - - - 

C  TEST  FOR  THE  OUTER  EDGE  OF  BOUNDARY  LAYER 

C  -  -  - - - - - - - - - - - - - - - 

C  COMPARE  THE  DIFFERENCE  OF  FLOW  PROPERTIES  BETWEEN  ThC  LAST 

C  TWO  CONSECUTIVE  GRID  POINTS  WITH  THE  SPECIFIED  TOLERANCE 

C  —  -  — - - - - - - - - - - - - - - 

ITEST=0 

140  ERROR=DABS (U(N)HJ(Nl) ) 

IF  (ERROR  LT  EPS)  GO  TO  150 


•(  '» 


»v* 

i 


C  -  '  - - ________________________ - - 

C  ADD  A  NEW  GRID  POINT  AT  THE  NEW  GRID  LINE  (M+l) 

C  -  -  - - - - - - - - -  - - - - -  -  -  - - - 

I TEST= ITEST+1 
U ( N+l ) =U<  N  > 

V ( N+ 1 ) =V<  N  > 

T(N+1)=T(N> 

RO<N+l )=RO<N> 

UP (N+l )=UP (N) 

VP  < N+l  )  =VP  <  N ) 

TP  (N+l  )=JP <N> 

ROP(N+l )=ROP(N) 

C  -  -  - - -  -  -  -  -  -  -  -  ________  ___________ 

C  CALCULATE  THE  MATRIX  COEFFICIENTS  AT  THE  NEW  GRID  POINT 

c_.  _________________________ - - - 

IF  (ITEST  GT  1)  GO  TO  160 
CALL  TH0M1  <N,  A,  B-  C,  D> 

CALL  THOMS  <N,  U,  A,  B.  D> 

GO  TO  170 
160  MU  1  =  1  D+00 
DMUI=0  D+OO 
DEUI=0  D+00 
DDEUI =0  D+00 
ROMD=ROP(N  >*DMDB<N  > 

UDI=ROMD*UP  <N) 

CR (N)=CITA*ROMD 
CRB (N ) =C I T  A1*R0MD 
CALL  V  ALUS  <  KE .  KE2 ) 

CALL  COEP 1  <N,  U.  RUI,  RVI ,  MUI ,  DMUI ,  DEUI,  DDEUI ,  UDI,  CR,  CRB,  A,  B,  C,  D) 
CALL  TH0M1  < N,  A.  B,  C,  D) 

CALL  TH0M2<N,  U.  A,  B.  D) 

170  N=N+1 

N 1 =N1 + 1 
N2=N2+1 
RU(N)=RU(N1 ) 

RUB (N ) =RUE (N1 ) 

RVB (N) =RVB (N1 ) 

DMDB ( N ) =DMDE ( N 1 ) 

DMNB ( N  >  =DMNB ( N 1 ) 

UD  < N) =UD ( N 1  ) 

VD (N)=VD(N1 ) 

TD ( N) =TD( N 1 ) 

IF  (N  EG  100)  GO  TO  999 
GO  TO  140 

150  IF  (ITEST  EQ  0)  GO  TO  200 
DO  180  1=1.  N2 
J=N-I 

1 80  CALL  THOMS  <  J,  U,  A,  B ,  D) 


D.6 


n  n  o  pj  onoooopjoonnnn 


SOLVE  THE  NEW  GAS  TEMPERATURE  T  AT  THE  NEXT  GRID 
LINE  ( M+l )  USING  SIX-POINT  IMPLICIT  SCHEME 


SET  THE  BOUNDARY  CONDITIONS  AT  THE  LINE  (M+l) 


OG  T ( 1 ) =TW 

T ( N  >  =  1  .  OD+OO 


CALCULATE  THE  INITIAL  VALUE  OF  SOME  COEFFICIENTS 


CALL  VALU1 (Cl  1 ,  CI2> CI3, CI4,  DY»  DY1 > 


CALCULATE  THE  COEFFICIENT  MATRIX  ELEMENTS 


DO  210  1=2,  N1 

CALL  PARA1  (  I.  U,  V,  RO,  RUI,  RVI  ) 

RUI=CITA*RUI+CITA1*RUB  <  I) 

CALL  PARA3  <  I,  KE2,  KE3,  U,  DEUI  > 

CALL  P  ARA3  (I,  KE2,  KE3,  T,  DETI) 

DY1=DY1*KE 
DU I=DEUI/D Y 1 
DT I=DETI/DY 1 

CALL  PARA3(  I,  KE2,  KE3,  RU,  DEUI ) 

DUI=CITA*DUI+C ITA1*DEUI/DY1 
CALL  PARA4(  I,  KE,  KE1,  T.  DDETI  > 

CALL  PARA2  (  I,  T,  MUI ,  DMUI,  MDR,  MNU) 

UP  U=C I T A*  <  UP ( I  >-U< I )  >+CITAl*(UP(I >-RU< I ) ) 

UD I=ROP < I ) *UPU*DMDB ( I > 

EDM=EC  *DU I *DU I *MU I *DX 
ERD=EC*UDI*UPU 
ROMN=ROP( I >*DMNB( I ) 

TPT=R0MN*TP(I > 

TD I =EDM+ER  D+TP  T 
CR (  I ) =C ITA+RQMN 
CRB ( I >  =CITA1*R0MN 
CALL  VALU2 (KE,  KE2 ) 

MU  I =MU I /PR 
DMUI=DMUI*DTI /PR 

10  CALL  C0EF1  (  I,  T,  RUI ,  RVI ,  MUI,  DMUI,  DETI,  DDETI ,  TDI .  CR,  CRB,  A,  B,  C,  D) 
CALL  C0EF2  (  T ,  A,  B,  C,  D) 


SOLVE  THE  GAS  TEMPERATURE  PROFILE  BY  THOMAS  ALGORITHM 


DO  220  1=2, N1 
;  220  CALL  TH0M1  (  I,  A,  B,  C,  D) 

|  DO  230  1=1,  N2 

'  J=N-I 

CALL  TH0M2<  J,  T,  A,  B,  D) 


230 


<o  o 


C  CALCULATE  THE  GAS  DENSITY  PROFILE  AT  THE  GRID  LINE  <M+1 ) 

C  -  -  -  --  --  --  --  --  --  --  --  --  --  --  --  --  --  - 
DO  240  1=1 ,  N 
240  R0(  I  )=1.  D+OO/T  ( I ) 

C  -  -  -  -  -  -  -  -  -  -  -  -  -  - - -  - - - 

C  SOLVE  THE  GAS  NORMAL  VELOi  ITY  PROFILE  AT  THE  GRID  LINE  (M+l) 

___________________________ 

DO  310  1=1  , N 
310  RU< I)=RO< I )*U(  I  ) 

DO  320  1=2, N 

B  <  I  >=RU<  I  > -RUE < I >  +  RU< 1 -1 > -RUB (  1-1 > 

32:  C  <  I  )=R  VD  I  )  -R'.'B  ( I  -  1  ) 

CF=DY/ ( C I 3*KE > 

D  <  1  )=0  D+00 
DO  330  1=2  N 
CF=CF*KE 

22:  Dl  I  )=D (  1-1  >-CI TA2«C ( I ) -CF*B ( I ) 

DO  340  1=1 , N 
340  V ( I )=D ( I ) *  T ( I ) 

r  -  ______________________________  _ 

C  SOLVE  THE  NEW  PARTICLE  PARAMETERS  ( UP,  VP,  TP,  AND  ROP  )  AT  THE 

C  GRID  LINE  <M+1>  USING  FOUR-  OR  SIX-POINT  IMPLICIT  SCHEME 

c_____________________  ___________ 

C  CALCULATE  SOME  COEFFICIENTS  OF  THE  FINITE  DIFFERENCE  EQUATION 

c_._ - - - - - - - - - - - - - - - 

DO  410  1=1 , N 
UDI=UDi  I  ) 

VD I =VD ( I ) 

TD I =TD  <  I  ) 

CALL  PARA2 (IT, MUI , DMU I .  MDR.  MNU ) 

CALL  PARA5  (  I,  DX,  MDR,  MNU,  DMD,  DMN,  U,  V,  T,  UD,  VD,  TD  ) 

UD< I )=CITA*UD( I )+C ITA1*UDI 
VD ( I )=CITA*VD< I )+C ITA1*VDI 
TD ( I)=CITA*TD< I >+CITAl*TDI 
DMD ( I ) =C I T  A*DMD ( I ) 

DMN (I ) =CITA*DMN< I ) 

DMDB ( I  )=C I TAi*DMBB  < I ) 

410  DMNB( 1 )=CITA1*DMNB < I ) 

c_._________________________ - - 

C  STORE  THE  VALUES  OF  PARTICLE  VELOCITIES  AT  THE  PREVIOUS  LINE 

C  -  -  - - __________ - __________ - - - 

DO  420  1=1 , N 
RUB  <  I  >=UP<  I  ) 

420  RVB  ( 1  )  =VP  (  I  ) 

C  _  _  - - - - - - -  - - -  - - - - - - - - 

SOLVE  THE  PARTICLE  TANGENTIAL  VELOCITY  UP  AT  THE  LINE  <M*-1> 


IF  (X.LT  XCRI)  GO  TO  530 
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c 

c 


SOLVE  THE  DIFFERENCE  EQUATIONS  USING  SIX-POINT  SCHEME 


K 


C  SET  THE  BOUNDARY  CONDITIONS  AT  THE  LINE  <M+1> 

C  -  -  _______________________ 

UP  <  1  )=0  D+OO 
UP  (N)  =  l  D+OO 


C  CALCULATE  THE  INITIAL  VALVE  OF  SOME  COEFFICIENTS 


CALL  VALU1  <  C  T 1 »  C 1 2  CIO  -  CIA,  DY.  DY1 ) 


CALCULATE  THE  COEFFICIENT  MATRIX  ELEMENTS 


DO  510  1=2. N1 
RU I =UP i 1 ) 

R  V  I  =VP <  I  ) 

CALL  PARA3  <  I,  KE2,  KE3>  UP,  DEUPI  > 

DDEUP I =0  D+OO 
MU  1=0  D+OO 
DMUI=0  D+OO 
UPDI=UD( I ) 

CALL  VALU2 (KE, KE2 ) 

510  CALL  C0EF1  (  I.  UP,  RU  I ,  RV  I ,  MU  I ,  DMUI ,  DEUPI,  DDEUP  I,  UPD I-  DUD.  DMDB, 
1  A,  B,  C,  D  > 

CALL  C0EF2  ( UP  >  A,  B,  C ,  D  ) 

DO  520  1=2, N1 

520  CALL  TH0M1  (  I,  A,  B,  C,  D) 

DO  525  1=1,  N2 
J=N-I 

525  CALL  TH0M2  (U.  UP.  A,  E,  D) 

GO  TO  560 

C  -  -  -  --  --  --  - -  -  - 

C  SOLVE  THE  DIFFERENCE  EQUATIONS  USING  FOUR-POINT  SCHEME 

c  -  -  -  —  -  -  -  -  -  -  -  -  -  -  -  - - - - 

530  DX 1=C I 1*KE/DY 
DX2=CI2*KE/DY 
DO  540  1=1,  N1 
DX 1 =DX 1 /KE 
DX2=DX2/KE 
A ( I )=DX1*VP < I ) 

D( I)=UP(I ) -A< I ) +DMD ( I > 

DX  VP=D  X2*VP ( I ) 

540  C< I >=-DXVP*UP< I  +  1 >  +  (UP  < I ) +DXVP -DMDB ( I ) >*UP ( I ) +  UD ( I ) 

UP (N)= 1  D+OO 
DO  550  1  =  1,  N1 
J-N-l 

550  UP  ( J)  =  (C(  J)-A(  J)*UP(U+1  ) )  /B  ( J) 


r>  o 


C  -  -  - - -  -  -  - - - - -  -  -  - - - - - - - 

C  CALCULATE  THE  VALUES  Of  UP  AT  THE  HALF  POINT  <M+©> 

560  DO  570  1=1.  N 

570  UD<  I  >=CITA*IF<  I  )+C  ITA1  *RUB  <  1  > 

C  -  -  ___________________ - _____  -  - 

C  SOLVE  THE  PARTICLE  NORMAL  VELOCITY  VP  AT  THE  LINE  (M+l) 

C  BY  MEANS  OF  THE  FOUR-POINT  SCHEME 

c_.____________________________ 

DX 1=CI 1+KE/DY 

DX2=CI2+KE/DY 

DO  610  1=2, N 

DX 1 =DX 1/KE 

DX2=DX2/KE 

B< I )  =  - D  X 1  *  VP  < I  ) 

A  (  I  >  =  UD'  I  )  -I-  (  I  )  -*-DMD  ( 1  ) 

DX VP=DX2*VP ( 1 ) 

6  10  C  <  I  )=<  UD<  I  )  -DX  VP-DMDB  <  I  )  )*VP(  I  >+DXVP*VP  <1-1  )+VD  ( I  ) 

VP  <  1  ) =0  D+00 
DO  620  1=2. N 

620  VP ( I )  =  (C (  I  )-B  <  I  )*VP  U-1))/A(I) 

c_________________________ - - - 

CALCULATE  THE  VALUES  OF  VP  AT  THE  HALF  POINT  <M+©) 


DO  630  1  =  1  , N 

630  VD  <I)=CITA*VF<  I )+C ITA1  *RVB  <  I  > 

C  -  -  - - ________ - ______ - - - - 

C  SOLVE  THE  PARTICLE  TEMPERATURE  TP  AT  THE  LINE  (M+l  > 

C  _  -  - - -  - - - - ____ - - - - - - 

IF  (X  LT  XCRI )  GO  TO  730 

C  -  -  -  -  -  -  -  -  -  -  -  -  -  - - - 

C  SOLVE  THE  DIFFERENCE  EQUATIONS  USING  SIX-POINT  SCHEME 

C  -  --  --  --  --  --  --  --  --  --  --  --  --  --  -- 

C  SET  THE  BOUNDARY  CONDITIONS  AT  THE  LINE  (M+l) 

c_._________________________ - 

TP  < 1 )=TW 
TP ( N ) = 1  D+00 

C  -  -  -  ____________________ 

C  CALCULATE  THE  INITIAL  VALUE  OF  SOME  COEFFICIENTS 


CALL  VALU1  (CI1.CI2.CI3, CIA, DY,  DY1) 


C  CALCULATE  THE  COEFFICIENT  MATRIX  ELEMENTS 

r  -  -  ___  ____  ______________ 

DO  710  1=2. N 1 
RU I =UD t  I ) 

RV I =VD < I ) 

CALL  PARA3  (  I.  KE2,  KE3.  TP,  DETPI  ) 

DDE  TP  I  =0  D+00 
MU  1=0  D+00 
DMUI-0  D+00 
TPDI=TD< I ) 
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CALL  VALU2  <KE,  KE2 ) 

■’JO  CALL  COE^l  (  I ,  TP,  RUI , RVI ,  MU1,  DMUI,  DETPI  .  DDE  TP  I .  TPDI .  DMN.  DMNB. 

A.  E,  C.  D> 

CALL  C OEr 2<Tf.A,B.  C,D) 

DO  720  1  =  2.  N1 

t:  CALL  THOMl  <  I,  A,  B.  C,  D> 

DO  725  1=1,  N2 
J=N-I 

’’O:  CALL  TH0M2<  J.  TP,  A.  B,  D> 

GO  TO  805 

c  _____________________ - ______ - - 

C  SOLVE  THE  DIFFERENCE  EQUATIONS  USING  FOUR-POINT  SCHEME 

C  -  -  -  --  --  --  ______  ----------------- 

■'30  DX  1=C  1  1  *KE  /Dy 
DX  2=C I 2*kE  /Dv 
DO  740  1=1 , N1 
DX 1 =DX 1 /KE 
DX  2=DX  2/KE 
A (  1  ) =D  a  1 *VD ( I ) 

B  <  I  )=UD(  I)-A(l  ) +DMN ( 1  > 

DXVP=DX2*VD<  1  ) 

741  C (  I  )  =  - D  X  VP  *  TF (  1  +  1 )  -MUD < I )  ■+  DXVP -DMNB i I )  > *TP ( I )+TD< I  ) 

TP ( N ) = 1  D+00 
DO  750  1=1,  N1 
J=N-I 

7  50  TP (J>= (C  <  J)-A( J>*TP (J+ 1  ))/B(J) 

r_.___________________________ - - 

C  SOLVE  THE  PARTICLE  DENSITY  ROP  AT  THE  LINE  (M+l) 

C  -  -  -  --  --  --  --  --  --  --  --  --  --  --  --  --  --  - 

905  IF  (X  LT  XCRI)  GO  TO  e30 

c_.___________________________ - - 

C  GET  THE  PARTICLE  DENSITY  BY  THE  ASSUMPTION  THAT  THERE  IS  NO 

C  DEPOSITION  OF  PARTICLES  ON  THE  SURFACE  OF  THE  PLATE 


DO  820  1=1,  N 
820  ROP< I ) =BETA*RO ( I ) 

GO  TO  900 

c_.___________________________ 

C  SOLVE  THE  DIFFERENCE  EQUATIONS  USING  FOUR-POINT  SCHEME 

r  -  -  ___________________________ 

63I-  DX  1=C  I  1*KE /DY 
DX2=C12*KE/DY 
DO  840  1=1, N1 
DX 1=DX 1/KE 
DX  2=DX  2/KE 
A ( I )=DX1*VD< I ) 

B ( I )=C1TA3*UP( I ) -C I TA4  *RUB ( I ) 

1  +C ITA*DX 1  * (VP ( I+l ) -2.  D+00*VP< I ) ) 

2  +C ITA*DX2* (RVB  <1+1 )-2.  D+00*RVB(I>  > 

840  C< I )=-DX2* VD( I ) *ROP ( 1+ 1 )+ ( C ITA4*UP ( I ) +  C ITA5*RUB ( I ) 

1  -CITA1*DX 1#<  VP  <1+1 )—2.  D+00*VP(I>> 

2  -CITA1*DX2*<RVB<  I  +  I  >-2.  D+00*RVB<I >  > )#ROP( I ) 


D.ll 


ROP (N ) =BET  A 
DO  850  1=1.  N1 
J=N— I 

85 j  ROP(J)=<C( U)-A< J)*ROP< J+l ) )/B( J) 

900  CONTINUE 

c_______ - - - - - -  - - - 

C  CALCULATE  THE  CHARACTERISTIC  QUANTITIES  OF  BOUNDARY  LAYER 

C  FLOW  SHEAR  STRESS.  HEAT  TRANSFER  AND  DISPLACEMENT  THICKNESS 

r  -  -  ___  _  _  _  _  __  __  ____  ____  _____  __  _____ 

C  GET  THE  SHEA9  STRESS  AT  THE  WALL 


CMDY=CK1*MUW/DY 

SHEAR= ( U( 2  >  — IM  1 ) )- <U(3)-U( 1 ) )/C K2+ (U(4 )-U( 1 ) ) /CK3 
SHEAR = SHEA 9 *CMDY 
r--—— 

C  GET  THE  HEAT  TRANSFER  AT  THE  WALL 

r  -  -  _______________  ________  _  _  _____ 

CMDY=CMDY/ <PR*EC> 

HEAT=( T(2)-T( 1 ) )- ( T <3 ) -Y< 1 ) ) /CK2+ ( T < 4 ) -T ( 1 ) >/CK3 
HE AT=-HEAT  *CMDY 

C  _  -  - - - - - __________ - -  -  -  - - - - 

C  GET  THE  DISPLACEMENT  THICKNESS  USING  THREE -POINT  DIFFERENCE 

C  FORMULA 

c_.  ___________________________ - - 

DYl=DY/(fc  D+00*KE) 

DO  910  1  =  1 .  N 
910  RU(I)=1  D+OO-RUII) 

RU (N+l ) =0  D+00 
TH ICK=0  D+00 
DO  920  1  =2.  N 
DY1=DY1*KE 

SUMI=CK4*RU< 1-1 )+CK5*RU< I ) -CK6*RU< 1+1 ) 

SUMI=SUMI*DY1 
9?0  TH I CK  =  THI C  k  +SUMI 

t  -  -  - - -  -  -  - - - ______ - - - - - 

C  OUTPUT  THE  COMPUTATION  RESULTS  AT  THE  GRID  LINE  (M+l) 

c_._____ - ____________________ - - 

950  WRITE<6,  1 1  10)  M.  N.  X,  SHEAR.  HEAT,  THI  CK 
WR  I TE  (  8 ,  1100)  <U(I  > .  V<  I),  T  (I).  RO(I  ), 

1  UP  (  I )  >  VP  <  I  )  ,  TP  <  I  ) »  R  OP  (I  ) .  I  =  1 .  N ) 

IF  (X-XMAX )  10,  10,  999 
999  STOP 
END 


D.  12 


o  r»  o 


SUBROUTINE  PROGRAMS 


C 

C 

c 


C  -  -  -  —  _______________________ 

SUBROUTINE  PAR A1  <  1 .  WJ>  WV>  WR ,  RU I »  RV I  > 

C  CALCULATE  SOME  PRODUCTS  USED  WHEN  SOLVING  THE  FINITE 

C  DIFFERENCE  EQUATIONS 

r  _  _  ^  _  _  _  _  _  _  _  _  _  _  _  _  — .  —  _  . 

IMPLICIT  REAL*B(A-H,  O-Z  > 

DIMENSION  WU(  1  ),  WV  <  1  ),  WR  <  1  ) 

C 

RUI=WR  <  I  >*W U(  I  ) 

RVI=WP  <  I  >*WV<  I  ) 

r 

RETURN 

END 


SUBROUTINE  PAR  A2 I,  T.  MU1.  DMUI .  MDR ,  MNU  ) 


C  CALCULATE  SOME  PARAMETERS  USED  WHEN  SOLVING  THE  FINITE 

C  DIF cERENCE  EQUATIONS 

c_.  - 

IMPLICIT  REAL*8(A-H. O-Z ) 

RE  AL*8  MU  I  ,  MDR  ..  MNU ,  NU 
DIMENSION  T  < 1  ) 

COMMON  /G/  PR,  WN,  S,  REYP,  EC 

C 

MU  I  =T  <  I  )  *■*  WN 
DMUI=WN*MUI/7 ( I ) 

DR  =  1  D+OO 
NU=2  D+OO 
MDR=MUI*DP 

MNU=MUI*NU/ <3  D+00*PR) 

C 

RETL^N 


SUBROUTINE  PArA3(  I  .  KE2,  KE3,  W.  DEWI  ) 
r  -  .  _  _  _  _  _  __  __  _  _  _  _  _  _  ______________ 

C  CALCULATE  THE  FIRST  ORDER  DIFFERENCE  USED  WHEN  SOLVING  THE 

C  FINITE  DIFFERENCE  EQUUAT I 0N5 

C  -  - - ______________ - - 

IMPLICIT  REAL*8(A-H,  O-Z  ) 

RE AL*8  HE2, KE3 
DIMENSION  W< 1 ) 

C 

DEWI»W< I  +  l  )+KE3*W( I )-KE2*W< 1-1  ) 

C 

RETURN 

END 
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C  -  -  - - - - - - - - - 

SUBROUTINE  PARA4< I , KE.  KE1 .  W, DDEWI ) 


CALCULATE  THE  SECOND  ORDER  DIFFERENCE  USED  WHEN  SOLVING  THE 
FINITE  DIFFERENCE  EQUUATIONS 


IMPLICIT  REAL*B(A-H,  0-Z  ) 

RE AL*8  KE,  KE1 
DIMENSION  W ( 1 ) 

DDEWI  =W  <  I-M  )-KEl*W  (  I)+KE*W<  1-1  ) 

RETURN 

END 


SUBROUTINE  PARA5C  I ,  DX,  MDR.  MNU.  DMD,  DMN,  U,  V,  T,  UD.  VD.  TD) 


CALCULATE  SOME  PARAMETERS  RELATED  TO  THE  INTERACTION  TERMS 
BETWEEN  GES  AND  PARTICLES 


IMPLICIT  REAL*8<A-H,  O-Z  ) 

REAL*8  MDR , MNU 

DIMENSION  DMD  (  1  ),  DMN<  1  > ,  U<  1  > ,  V  <  1  ) ,  T  ( 1  )  ,  UD  (  1  ) .  VD  ( 1 )  .  TD  (  1  ) 

DMDCI )=DX*MDR 
DMN ( I >=DX*MNU 
UD ( I ) =U  < I ) *DMD  < I ) 

VD  < I ) = V ( I ) *DMD  < I ) 

TD ( I )=T( I ) *DMN ( I ) 

RETURN 

END 


SUBROUTINE  VALU1  (C  I  1.  C  12.  C  13,  C  14,  DY,  DY  1  ) 


GIVE  THE  INITIAL  VALUES  OF  SOME  COEFFICIENTS 

C  -  -  - - - - - - - - - - - -  ■ 

IMPLICIT  REAL*8(A-H.  O-Z  ) 

RE AL*8  KE.  KE1 

COMMON  /HI/  KE , KE 1 , KE2 , KE3 
COMMON  /H2 /  CA.CB.CC.CD 
C 

DY 1=KE 1 *DY 
DY2=DY1*DY/KE 
CA=C 1 1 /DY 1 
CB«=C12/DY1 
CC-CI3/DY2 
CD«CI4/DY2 
C 

RETURN 

END 
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c-------------- 

SUDROUT INE  VALU2 <  KE, KE 2 ) 


C  CALCULATE  THE  VALUES  OF  SOME  COEFFICIENTS  AT  THE  GRID  POINT 

IMPLICIT  REAL«8<A~H,  0-2  > 

REAL*B  KE  -  KE2 

COMMON  /HZ/  CA.CB.CC.CB 

r 

C A=CA/KE 
CB=CB/KE 
CC=CC / KE2 
CD=CD/KE2 

r 

RETURN 

END 

c_________________________ - _____ 

SUBROUTINE  COEFi  (  I  ,  W.  RUI.  R  VI,  MUI,  DMUI .  DEW  I ,  DDEWI.  WDI,  CR,  CRB, 
♦  A,  B.  C,  D> 


CALCULATE  THE  COEFFICIENT  MATRIX  ELEMENTS  A,  B,  C,  AND  D  FOR 
C  THE  CRANK-NICOLSON  SCHEME 

IMPLICIT  REAL*8(A-H,  0-Z  ) 

RE AL*B  KE.  KE1 .  KE2,  KE3,  MUI 

DIMENSION  W  <  1  )  .  CR  <  1).CRB(1  )<A(  1  ) ,  B  (  1  > ,  C  ( 1  ) >  D( 1  ) 

COMMON  /HI/  KE. KE1 , KE2, KE3 
COMMON  /H2  /  CA.CB.CC.CD 
C 

RVT=RV  I  -DMUI 
CR  VT=C  A*R V  T 
CMU=CC  *MU I 
A ( I )=C RVT-CMU 

B  (  I  >=RUI+HE3*CRVT-*-kEI*CMU+CR(  I  > 

C (  I  )  =  -KE2* CRUT -KE * CMj 

DEWI=DEWI*CB 

DDEWI  =  DDEW I *CE 

D (  I  )=( RUI-CRB  <  I  )  )*W( I ) -RVT*DEWI +MU I *DDEWI  +WDI 
C 

RETURN 

END 
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SUBROUTINE  C0EF2<  W,  A,  B  ,  C,  D  > 


I 

SPECIFY  THE  COEFFICIENT  MATRIX  ELEMENTS  A(1),B(1).C(1).D(1) 


IMPLICIT  REAL*8(A-H,  0- Z  > 

DIMENSION  W ( 1  ) <  A< 1  ) .  B  (  1  ) «  C  <  1 ) .  D(l) 

A  (  1  )  =  0  D+00 
B  (  1  )  =  1  D+00 
C ( 1 )=0  D+00 
D ( 1 >=W< 1 ) 

RETURN 

END 


SUBROUTINE  THOM  1  (  1  ,  A,  B  ,  C.  D  ) 


ESTABLISH  UPPER  TRIANGULAR  MATRIX  (FORWARD  ELIMINATION)  FOR 
SOLVING  THE  TRI-DIAGONAL  SYSTEM  FOLLOWING  THE  THOMAS  ALGORITHM 


IMPLICIT  REAL*8<A-H, 0-2  ) 

D I  MENS  ION  A(1)>B(1)>C(1)>D(1) 

BC=C(  I  )  /B  (  I  -1  ) 

B (  I )=B ( I )-BC*A (  1-1  ) 

D ( I )=D ( I >-BC«D< 1-1  ) 

RETURN 

END 


SUBROUTINE  THGM2<  I ,  W.  A ,  B.  D  ) 


GET  SOLUTION  (BACK  SUBSTITUTION)  FOLLOWING  THOMAS  ALGORITHM 


IMPLICIT  REAL*8(A-H,  0- Z  > 
DIMENSION  W(1).A(1).B(1)/D(1) 

C 

W  <  I  )  =  (  D  ( I  )  -  A(  I  )  *W  <  I  + 1  )  )  /B  <  I  > 

C 

RETURN 

END 
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APPENDIX  E 


AN  ADDITIONAL  DISCUSSION  REGARDING  THE  ASSUMPTION  OF  THE  PARTICLE-DENSITY 

PROFILE  AFTER  THE  CRITICAL  POINT 


In  the  foregoing  analysis,  it  was  pointed  out  that  at  the  critical 

point,  the  particle  velocity  at  the  wall  vanishes.  Physically,  upw  =  0 
means  that  the  particles  will  stop  and  stay  at  the  wall.  Consequently,  a 
lot  of  particles  would  gradually  accumulate  on  the  surface  and  it  would 

result  in  a  great  increase  of  the  particle  density  at  the  wall  unless  a 

diffusion  mechanism  is  the  predominant  process.  Therefore,  contrary  to  a 
diffusion-controlled  process  (as  assumed  in  the  previous  analysis),  the 
accumulation  of  particles  on  the  surface  represents  another  extreme  limiting 
case  for  a  gas-particle  system.  It  is  not  the  purpose  of  this  discussion  to 
obtain  information  concerning  particle  deposition  on  the  flat  plate, 

although  such  data  can  be  used  in  some  practical  problems  (for  example, 
retardation,  accumulation  and  impingement  of  particles  on  solid  surfaces 
i  ave  considerable  effects  on  the  erosion  of  the  surfaces).  The  aim  here  is 
to  attempt  to  gain  more  insight  into  the  assumption  of  the  particle-density 
pr  ofile  when  x  >  xcr-j . 

In  the  following,  it  is  assumed  that  the  accumulation  of  the  particles 
at  the  wall  is  allowed  and  the  thickness  of  the  particle  accumulation  layer 
can  be  neglected  compared  with  the  boundary-layer  thickness.  The  particle 
density  will  eventually  become  very  large  at  the  wall.  So  it  is  reasonable 
to  specify  the  reciprocal  of  particle  density  at  the  wall  to  be  zero.  Then 
an  additional  boundary  condition  is  obtained  and  the  six-point  scheme  can  be 
applied  to  the  continuity  equation  of  the  particle  phase.  Using  the 
reciprocal  of  the  density  as  a  new  dependent  variable,  the  continuity 
equation  (2.16)  can  be  written  as 


UP  I -  (-)  *  v„  JL  (L)  - 

p  dx  pp  p  ay  Pp 


(E-l) 


With  the  quotient  expressions  for  the  six-point  scheme  (3.1)-(3.4),  its 
finite-difference  form  is 
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where 


/D  9Avn  ,  +  (1-8)Avd 

pm+8,n  pm+l,n  pm,n 


The  boundary  conditions  at  the  wall  and  the  outer  edge  of  the  boundary-layer 
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By  solving  the  difference  equation  (E-2)  with  the  boundary  conditions  (E-7) 
and  (E-8),  a  numerical  solution  is  obtained  for  the  same  freestream  and  wall 
conditions.  In  Figs.  E-l  to  E-4,  these  results  are  compared  with  those 
results  obtained  based  on  the  assumption  of  pp  =  Pp  for  x  >  xcn-.  It  is 
interesting  to  note  that  the  resulting  boundary-layer  characteristics  (i.e., 
shear  stress,  wall  heat-transfer  and  displacement  thickness)  are  almost  the 
same  for  these  two  different  methods,  despite  the  fact  that  there  are  large 
differences  between  the  particle  density  profiles,  especially  near  the  wall. 
It  provides  the  evidence  that  the  particle  density  has  a  very  small  effect 
on  dusty-gas  boundary-layer  flows  in  the  quasi-equi 1 ibrium  and 
near-equi 1 ibrium  regions.  As  mentioned  before,  this  approximate  treatment 
of  the  particle  density  (that  is,  for  x  >  xcr-j,  Pp  =  Pp  is  assumed  across 
the  whole  boundary  layer)  results  in  satisfactory  accuracy  on  the 
boundary-layer  flow-profiles  and  characteristic  quantities  except  the 
particle  density  itself.  It  is  due  to  the  fact  that  the  particle  and  gas 
phases  are  already  in  a  near-  or  quasi-equilibrium  state  after  the  critical 
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